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ABSTRACT 

Pair annihilation of heavy stable particle that occurs in the early universe is 
investigated, and quantum kinetic equation for the momentum distribution of the 
annihilating particle is derived, using the influence functional method. A bosonic 

ON ■ 

field theory model is used to describe the pair annihilation in the presence of decay 
product particles making up a thermal environment. A crossing symmetric Hartree 
approximation that determines self-consistently the equilibrium distribution is de- 
veloped for an otherwise intractable theory. The time evolution equation and its 

Markovian approximation is derived, to give a generalized Boltzmann equation in- 

■ 

eluding off-shell effects. The narrow width approximation to an energy integral in 

C$ ' 

this equation gives the usual Boltzmann equation in a thermal bath of light particles. 
The off-shell effect is a correction to the Boltzmann equation at high temperatures, 
but is dominant at low temperatures. The effect changes the equilibrium distribu- 
tion from the familiar l/( e Wfc//T — 1 ) to a modified one given by a Gibbs formula. 
Integrated over momenta, the particle number density becomes roughly of order 



(coupling) x yT/M ■ T 3 at low temperatures for the S-wave annihilation. The relic 
mass density in the present universe is insensitive to the coupling strength in a large 
range of the mass and the coupling parameters, and scales with the WIMP mass as 
w 6 x 10 4 eV cm~ 3 (M / GeV) 4 ^ 3 . The bound from the closure density gives an upper 
WIMP mass bound roughly of order 1 GeV in the present model. 
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1 Introduction 



The Boltzmann equation for the particle distribution function is described in 
terms of S-matrix elements defined on the mass shell, such as the cross section and the 
decay rate. This integro-differential equation is intuitively appealing since relation 
to the classical theory is evident. In the quantum mechanical context its foundation 
and its possible generalization have extensively been discussed in the past. We may 
only mention two approaches; the real-time thermal Green's function method |1[ and 
the closed time path method |J . Both use as the fundamental quantity the Green's 
function that contains information off the mass shell. It is generally believed that 
some sort of coarse graining or reduction of detailed quantum mechanical information 
is needed to derive a quantum version of the Boltzmann equation. There is however 
no unique generalization of the Boltzmann equation known to us. Perhaps correction 
to the Bolzmann equation differs depending on a physical situation one has in mind. 

We investigate in the present work the non-equilibrium process of pair annihila- 
tion of some heavy particles in the presence of a thermal bath of lighter particles. 
Time evolution of the occupation number for the heavy particle in a particular mo- 
mentum mode is derived and analyzed from a new perspective. 

The problem is of great interest in cosmology, and occurs for pair annihilation of 
stable particles such as the anti-nucleon and the positron. For instance, the electron- 
positron annihilation has a great impact on nucleosynthesis that takes place imme- 
diately after this annihilation. The annihilation rate has been calculated using a 
thermally averaged Boltzmann equation, but estimate of the off-shell effect to the 
annihilation process has never been worked out (see, however, || for a limited calcu- 
lation of finite temperature effects). This list of interesting candidates for the cosmic 
pair annihilation should be expanded to include WIMP (weakly interacting massive 
particle) or LSP (lightest supersymmetric particle) [[| that has been hypothesized, 
with a good motivation, as a solution to the dark matter problem in cosmology. 

We would like to obtain a suitable and useful generalization of the Boltzmann 
equation in these circumstances. The physical situation we have in mind is charac- 
terized by the presence of annihilating heavy particles (e^ in the case of electron- 
positron annihilation) interacting among themselves (e^ scattering and pair annihi- 
lation into two photons) and with lighter particles (scattering with photons). Al- 
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though we work out a particular toy model for this, our method is general enough 
for extension to other cases, and moreover our result is simple enough for practical 
use. 

We believe that for derivation of the generalized Boltzmann equation it is impor- 
tant to use a general and flexible framework allowing for a consistent approximation 
scheme. It is often practically difficult to estimate correction to the Boltzmann equa- 
tion in the Green's function or the operator method. Our experience in a similar, 
but a much simpler problem of the unstable particle decay in thermal medium [Q 
— f?J suggests that the influence functional method invented by Feynman and Ver- 
non |H may be useful to the present problem. The quantum mechanical problem 
of the excited state decay 0, analogous to, but simpler than, the unstable particle 
decay in field theory, is completely solvable at the operator level , but at the same 
time the influence functional method gives an identical result || to the operator 
approach. In the present work we develop a new Hartree type of approximation of 
the annihilation-scattering problem within the influence functional approach. We 
also utilize the operator identity to facilitate our analysis. 

In order to quantitatively discuss the role of scattering in the annihilation process, 
it is important to simultaneously deal with the two processes, the annihilation and 
the scattering. It is thus best to employ a fully relativistic field theory that respects 
the crossing symmetry. For simplicity we take throughout this paper an interaction 
Lagrangian density of the form, <p 2 x 2 , where ip is a heavy bosonic field and x a lighter 
bosonic field. We assume that the mass of ip 3> the mass of x- The annihilation 
channel (pip —>■ XX is related to the scattering channel (px — ► <PX by a crossing 
symmetry. The inverse processes to these and the other 1 to 3 processes, ip <-> <pxx 
and x WXi that may also occur in thermal medium of x particles are treated here 
symmetrically. It is important to recall that a finite time behavior of the quantum 
system in thermal medium allows the process such as ip <-> tpxx > even if it is 
kinematically forbidden for the on-shell S-matrix element. 

We first derive an integral equation that self-consistently determines the equi- 
librium distribution function f(k) within the Hartree approximation. Derivation 
of this equation is based on the assumption that time variation of the distribution 
function proceeds more slowly than indivisual microscopic reactions occur. Under 
this assumption the result 0, |J of the completely solvable model may be used, and 
the large time (3> relaxation time ) limit can be taken. The resultant equilibrium 
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distribution function deviates from the ideal gas form, l/(e /3a;fc — 1), but may be 
understood by the Gibbs formula e~ ^ tot with H tot the total Hamiltonian including 
interaction between the ip system and the x environment. 

We then derive a time evolution equation for the distribution function. This 
equation contains the initial memory term, hence is non-Markovian, as any exact 
treatment of the quantum mechanical behavior would demand. We next devise a 
useful Markovian approximation to the Hartree model and examine this approxima- 
tion in detail. 

The on-shell Boltzmann equation arises as a result of the narrow resonance ap- 
proximation for a Breit-Wigner type of energy integral in our Hartree-Markovian 
model. The Boltzmann approximation is excellent in many practical cases, but it 
fails when the main contributing part to the energy integral includes the region off 
the resonance pole at uo the energy. Thus, the Boltzmann equation is modified 
significantly for the environment temperature <C the ip mass M. We present a com- 
plete kinetic equation for the momentum distribution function, which may be used 
at low temeratures. 

Integration over the momentum of the distribution function gives a rate equa- 
tion for the number density which is of prime interest in the annihilation-scattering 
problem. One may naively expect that the scattering process conserves the particle 
number, thus scattering terms cancel in the Boltzmann equation. This is explicitly 
confirmed for the on-shell part of our quantum kinetic equation. But it is not clear 
whether a similar complete cancellation occurs for the off-shell quantity. This is 
because this conservation law is based on the commutability of the particle number 
with the scattering part of the effective Hamiltonian, along with the unitary evolution 
in the quantum system. Dissipation due to the thermal medium and its associated 
fluctuation however causes a non-unitary evolution for a part of the entire system. 
It is thus an interesting open question whether the scattering-related contribution 
remains for the off-shell part. In any event we find that the scattering and its inverse 
process do contribute, but with a very small rate. A large number density of order, 
O[10~ 3 ] A \J~-^ T 3 , is derived for the off-shell contribution of the inverse annihilation 
process, where A is the relevant y? 2 x 2 coupling. At low temperatures this becomes 
much larger than the usual one, ~ (MT/2ir) 3 / 2 e~ M l T ; namely the thermal number 
density of zero chemical potential. 

In cosmology the thermal environment gradually changes according to the adia- 
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batic law; the temperature oc the cosmic scale factor l/a(t). Along with the obvious 
change of the number density oc l/a 3 (t), the thermally averaged rate in the general- 
ized Boltzmann equation decreases as the temperature decreases. Thus, the decou- 
pling or the freeze-out of annihilation takes place roughly at the temperature when 
the thermal rate is equal to the Hubble rate ||10|| . What is left after the freeze-out 
is then the relic abundance of heavy stable particles. It is important to accurately 
estimate the relic abundance of WIMP or LSP for the dark matter problem. 

The off-shell effect considered in the present work gives a relic abundance of 
order Y = n^/n 7 ~ O[0.1] (^~ l ) 1 ^ 3 , with m p \ the Planck mass. Unlike the estimate 
based on the on-shell Boltzmann equation, this abundance does not suffer from the 
suppression factor e~ M l T s where Tf is the freeze-out temperature, usually much less 
than M. Moreover, the relic fraction Y is insensitive to the y? 2 x 2 coupling, if the 
coupling is not too small, that is if A > 9.3 x 1(T 5 (M/GeV) - 32 for 10~ 3 GeV < M < 
1 TeV. Translated to the mass density of the present universe, this gives an allowed 
WIMP mass range for the dark matter; M < 1 GeV. It is of considerable interest to 
work out the allowed parameter region for LSP in realistic supersymmetric theories. 

The rest of this paper is organized as follows. In Section 2 the influence func- 
tional method is briefly explained, taking our model of <p 2 x 2 interaction. A non- 
perturbative approximation of the Hartree type is then introduced in connection 
to the influence functional. The kernel function in the reduced Hartree model is 
determined in terms of the spectral function that itself contains the correlator to 
be determined. In Section 3 we explain how the concept of the slow variation of 
the particle distribution helps to derive a self-consistent equation that determines 
the equilibrium distribution function. This is achieved by using the result of the 
exactly solvable model of the excited state decay. In particular, the explicit form of 
the coincident time limit of two-body correlators is required for its derivation, and 
we calculate the correlator by using a new generating functional method within the 
influence functional formalism. 

In Section 4 the time evolution equation is derived for the occupation number. 
This equation is non-Markovian and the initial memory term is identified. A simple 
Markovian approximation becomes possible, again under the assumption that the 
time scale of the (p number density variation is larger than that of the microscopic 
reaction time. The Markovian equation thus derived makes its relation to the usual 
Boltzmann equation transpant. We then point out that this Markovian equation con- 
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tains the off-shell effect which becomes dominant at low temperatures. Cancellation 
of the scattering process in the on-shell contribution, when one integrates over mo- 
menta, becomes important to get a large off-shell contribution to the number density 
at low temperatures. In Section 5 the cosmological estimate of the relic abundance 
of heavy stable particles is given in the present model. The off-shell effect prolongs 
the freeze-out epoch, but due to a power law decrease (ex T 3 - 5 ) of the number density 
with temperature, the net effect gives a larger freeze-out density. Finally a simple 
estimate of the present mass density is given; p ~ 0[6 x 10 4 ] (M / GeVY^ 3 eV cm~ 3 
in a certain range of the parameter space (M , A) . 

In four Appendices we explain technical points relegated in the main text; (A) 
generating functional method in the influence functional framework, (B) renormal- 
ization of the off-shell distribution function, (C) some details for the unstable particle 
decay treated in the Hartree approximation, and (D) technical details for computing 
various integrals needed to obtain the off-shell distribution function. 

A short summary of our result stressing the estimate of relic abundance has been 



given in a previous note In the present longer paper we give all details stated 
there and present further new results not given there. 



2 Hartree approximation 

We consider a relativistic field theory model of the Lagrangian density given by 

C = C v + C x + C- m t j (2-1) 
+ C x = \ (d^) 2 -\m*^+ 1 - (d, X ) 2 -\™ 2 X\ (2.2) 

A„t = -^V- (2-3) 

It is later explained how to renormalize and introduce counter terms when it becomes 
necessary. The heavy particle (p can pair-annihilate into a % pah; (p<p — > XX with 
the dimensionless interaction strength A. This coupling A must be less than unity 
for our method to be useful. We assume that the lighter particle % m &kes up a 
thermal environment of temperature T = 1/13 in our unit of the Boltzmann constant 
ks = 1 ■ It is thus assumed that although not written explicitly, there is interaction 
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among x particles themselves or with some other light particles to maintain a thermal 
equilibrium. 

In the problem of our interest one focusses on a particular dynamical degree 
of freedom (the ip field in our case) and integrates out the environment part (the x 
field) altogether. In the influence functional method || this integration is carried out 
for the squared amplitude, namely for the probability function directly, by using the 
path integral technique. This way one has to deal with the conjugate field variable ip' 
along with (p, since the complex conjugated quantity is multiplied in the probability. 

Define first the influence functional T by integrating out the x fi e ld degree of 
freedom during a fixed time interval between t{ and tf] 



j Vx J Vx' exp i J dx (C x (x) - L y 

+ C mt (p(x) , X (x)) - C int (p\x) , x\x)) )] . (2.4) 



We convolute with this the initial and and the final density matrix of the x sys- 
tem. For the initial state it is assumed that the entire system is described by an 
uncorrelated product of the system and the environment density matrix, 

P { t x pS x) , pS x) =pp = e-^W/tre-^W . (2.5) 

Here pp is the density matrix for a thermal environment, written in the operator 
notation using the environment Hamiltonian Hq(x)- To avoid a possible confusion, 
we sometimes write the operator explicitly by the tilded letter such as x here. The 
density matrix py' for the ip system is arbitrary, except that it is assumed to commute 
with the (p Hamiltonian, 

[ PV °\H Q (<p)]=0. (2.6) 

We believe that this restriction on the <p state does not exclude many practical cases 
of interest. At the final time tf the x integration is performed taking the condition 
of non-observation for the environment, 

J d Xf J dx' f 5( X f-x'f) (•••), (2-7) 

with the understanding that the environment is totally unspecified at the time tf. 

The result of x integration is given by a Gaussian integral and the influence 
functional to order A 2 is of the form, 

Tt^p,tp'\ = 
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ex P[-T / dxdy (£ 2 (x)a R (x ,y)£ 2 (y) + i£ 2 (x) a^x ,y) X 2 {y) )] , (2.8) 

4 Jxo>yo 

X 2 (x) = <^ 2 (x) + ^ 2 (x) , = ¥> 2 (s) - ¥>' 2 (x) , (2.9) 

a(x , y) = a R (x , y) + iaj(x , y) = A 2 tr (T[x 2 (x)x%) P/3]) • (2.10) 

Note the presence of the time ordering, xq > yo, in the above formula. Unless a 
confusion occurs, we shall simplify the four dimensional integral such as / d 4 x (• • •) 
by writing / dx (•••)• We note that the kernel function a(x ,y) satisfies the time 
translation invariance, thus may be written as a(x — y), as can explicitly be shown 
by taking the complete set of H (x) eigenstates; 

tr (A(x ) B(y ) p p ) = e-^n-mxo-yo) ^\ n ) (n\B\i) . (2.11) 

n ,i 

An explicit form of the kernel function a(x) or its Fourier transform is given later. 

Higher order terms in A 2 are actually present in the exponent of the influence 
functional. These contribute either to many-body processes we are not interested 
in or to higher order terms in our process. Moreover, Feynman and Vernon proved 
that the above form of the influence functional satisfies the required fundamental 
property such as the causality and the unitarity. We may thus safely ignore these 
higher order contributions in the weak coupling limit. 

Another point is the tadpole contribution and how the kernel oti is modified due 
to the tadpole, which we discuss later. (See eq.( [2.60D for the modification.) 

The convolution with the system variable (p gives the reduced density matrix for 
the system at any given time tf] 

(2.12) 

from which one can deduce physical quantities for the (p system. Here S((p) is the 
action for the ip system obtained from the basic Lagrangian. 

It is often useful to introduce a notation for the correlator in the influence func- 
tional method; for x > y , 



(<p{x)<p(y)) = J d<p(x) J d<p'(x) J d<p(y) J d<p'(y) 
V V [ V V 'e^-^^, V '}^xMy)p^\^y),^( y )). (2.13) 



In this formula the functional integration is performed for ip during the time interval, 
xq > t > y , where tf > xq > yo > U. The path integral prior to the time yo gives the 
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reduced density matrix p^ at the time yo, evolved from p^ at the time t{. In general, 
the correlation function (<p(x)ip(y)) does not obey the time translation invariance, 
reflecting that the initial memory is never completely erased in the quantum system. 
A nice feature of this correlator formula is that the initial memory effect appears 
compactly via the reduced density matrix p( R >. One may generalize the concept 
of the expectation value (• ■ ■) to any multiple of local operators including also the 
conjugate <p'. 

The model thus specified is difficult to solve due to the appearance of the quartic 



term of (p in the influence functional ( |2.8|) . The situation is however simplified 



when one considers a mean field approximation. In the mean field or the Hartree 
approximation one replaces a product of multi-field operators by a product of two 
operators with several averaged two-body correlators. This approximation is good if 
one can ignore a higher order correlation than that of two-body. The Hartree model 
is expected to work well in a dilute system. The dilute system is defined by the low 
occupation number for each mode. It should not be confused by a possibly large 
value of the number density which is a mode-summed quantity. For the bose system 
we consider here this is a circumstance very far from the bose condensed state. The 
dilute approximation seems good in most cosmological application. 

The Hartree approximation we now introduce is a Gaussian truncation to the 
influence functional; we replace the original one by properly defining a new kernel 
function (3{x ,y) in the quadratic form; 

•T^bV] = exp[- / dxdy {£{x)p R (x,y)£(y) + i£(x)l3 I (x,y)X(y))], 

(2.14) 

X(x), {(x) = <p(x)±<f/(x). (2.15) 

As for the correlator, the new kernel $ in the truncated model does not satisfy the 
time translation invariance. The identification of the new kernel (3{x ,y) is made by 
comparing two point correlators to any arbitrary order A 2 ; 

{X(x)£(y)) , (X(x)X(y)) , (£(x)£(y)) • It turns out that the last correlator (£(x)£(y)) 
vanishes both for the original and the Hartree-approximated model, which can be 
regarded as a consistency check of our approach. Once the Hartree model is de- 
termined by the kernel (3(x,y), one can work out its consequences to all order of 
A. 
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Let us first introduce two types of propagator for (p field; 

G(x,y) = i{<p(x)<p(y))r = r A , (2.16) 

and % ((p(x)tp(y))jr = i . The first one is the full propagator taking into account the 
environment interaction. The second one, on the other hand, is an extension of the 
free propagator in perturbative field theory accommodated to a non-trivial ip state 
in our problem. Without the non-local interaction in the influence functional J 7 , one 
may use for i (ip(x)ip(y))^=i the complete set of eigenstates of the ip Hamiltonian 
H (ip), to show that this quantity is translationally invariant; introducing a new 
notation for the propagator, 

G (x -y) = i (ip(x)ip(y))f =1 
= i r dk e- ik0 ^-^ Yl S ( k o~ E n + Ei) (#(0 ,x)\n)(n\<p(0,$)\i). (2.17) 

J -°° n,i 

Note that pf^ is time independent. It is always useful to decompose this into inde- 
pendent momentum modes; for xq > 

-iG (x) = j j^j- (/o(^ + (l + 0))e- fa ) , (2.18) 

where ko = Uk = \Jk 2 + M 2 . One may use the harmonic oscillator basis (n = 
, 1 , 2 , • • •) in this plane wave decomposition to show 

/ (£)H<l|y5(£)|0>/£np?. (2.19) 

n 

In the usual language of the creation and the annihilation operators, 

(aia 2 )r =1 = = (a\al) r=1 , (a\k')a(k)) T=1 = f (k) 5(k - %') . (2.20) 

We compute correlators in two theories, using the two different forms of the 
influence functional, JF 4 and JF 2 , and identify two results. One then has an equality; 

J dx'dy'G (x-x')(3(x',y')G (y'-y) 
= -i J dx'dy'Go(x-x')a(x'-y')G(x',y r )G (y'-y). (2.21) 

This identity leads to a crucial relation of the two kernels; 

(3(x,y) = -ta(x-y)G(x,y). (2.22) 
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It is important in the self-consistent Hartree approximation that the full propagator 
G instead of Go appears in this equation, since the truncated model should contain 
the full <p propagator to be determined self-consistently. Equivalently, the full prop- 
agator has a self-consistency condition in the Hartree approximation; in the matrix 
form, 

G = G - i G (3G + • ■ ■ = G - G (aG)G + - • • . (2.23) 

The next task is to derive a useful relation between the spectral functions that 
appear in eq. (|2.22 ). The full propagator G has two spacetime arguments, for which 



we use the center of mass (CM) variable (x + y)/2 and the relative coordinate x — y. 
We observe that there is no dependence on the CM space coordinate (x + y)/2 due 
to the spatial homogeneity, hence Fourier-transform the relative coordinate to get 

- i G(x , y) = J |^ /(- k , ^±*) e*M . (2.24) 

The spectral / here may depend on the CM time. 

The occupation number f(k,t) that becomes important in subsequent discussion 
is defined from the coincident time limit of the Fourier transformed full propagator. 
We first define the spatial Fourier decomposition of the full propagator, 



G(x ,y ;k) = j d 3 (x-y) G(x ,x;y ,y) e ~ ik <^ (2.25) 

= i (ifi(k,x )ip(-k,y )) , 
(p(k,x )= [ d 3 xip(x ,x )e- ife . (2.26) 



Due to the hermiticity of the field ip, 

ft$,x )=<p(-Z,x Q ). (2.27) 

Using the Heisenberg equation for (p, one has for the occupation number 

- (tit \ , X \ - 1 / # t fe g g) d<f(k,x ) _ 2 _ ut r 

u k [f{k ,x ) + -J = -{ ^7 ^7 — + uJ k p J {k , xoMk , x ) ) (2.28) 



i 
4 



d 9 \2^, T\ 
- G(x ,yo;- k) 



dx dy ' 



~^(^(k,x )px 2 (k,x )}. (2.29) 



The first quantity in the right hand side is further related to / by 

1 1 r At, _ \ 

f{k,t) + - = — -^klf(k,t)- — (^(k,t)^(k jt )}. (2.30) 
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Here uJk is what we may call a reference energy to define the occupation number for 
the mode k. We take for the time being uJk = uJk, the energy of free particle, and 
later modify the definition of the occupation number f(k , t) slightly by allowing the 
temperature dependent (p mass in u k . 

We define the spectral weight for the two kernels in terms of the Fourier decom- 
position of the relative spacetime coordinate; for xq > yo 

= / (0 rdh^ r - {k) e ~ iHx ~ v} ' (2 - 31) 

ft* ,y) = f r(k , ^ ) e-*(-*) . (2.32) 

A strange looking factor 1/(1 — e _/3fco ) is inserted for the kernel a(x), since this factor 
compensates the non-analytic property of the real-time thermal Green's function 
a(x) defined in eq. (|2.10| ), relative to the analytic imaginary-time Green's function. 
The factor is well understood [|1^] and is also explained in our previous paper [13|. We 
shall shortly give an explicit formula of r x (k). The spectral weight r for the Hartree 
model is determined from eq.( |2.22| ); it is of the form of a convolution integral; 

/, x f d4k ' ?, >/ n 2rJk + k') 

r(M) = 7 (2 - 33) 

The four dimensional integral fl2.33|) is not very useful in practice. In the following 
we shall derive an important relation of the two spectral functions, r x and r, using 
the occupation number / and another quantity v; 



(27r) 3 2cj fc , I 1 - e-i 3 ^^) 

- i j_ # '(g, t) ^ ( g, f ) _ ^ t , - 

2 uj^ dt dt 
In order to derive this relation, we first note that the quantity of the form, 

/ ^f(-k ,x )F(k ), (2.36) 
J In 

appears in eq. ( |2.33|) . This is related to a coincident limit of the full propagator, 
dkn ~, 



j ^f(-ko,x )F(k ) 
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= lim F 



lim F I — 



i , d 
2^dx^ 



dy } J 2tt /l 



(2.37) 



For simplicity we omitted the spatial coordinate or its Fourier component. 

We expand in the Taylor series the derivative operation applied to 
— iG(x,y) = (ip(x)ip(y)) in this equation and evaluate derivatives using the Heisen- 
berg equation for ip. For instance, 
d d 

lim + ("a "a ) V&My) = <p(x,x )(p(y, xq) - ip(x , x )(p(y , x ) , (2.38) 

xo ->y+ ox dy 



lim 

xo^Vo VXq 



d d 

— ) 2 ip(x)ip(y) 



dy ' 



(V 2 - M 2 ) ip(x,xoMy, xo) + <p(x , rr )(V 2 - M 2 ) cp(y, x ) 
-2<p(x ,x )<p(y ,x ) - \((px 2 )(x,x )ip(y,x ) - \tp(x , x )((px 2 )(y , x ) . (2.39) 

These terms have different structure for the even and odd number of derivatives, 
and it is convenient to separately deal with the even and the odd powers of the 
F(u) expansion, with the exception of the first term F(0). One thus arrives at the 
following operator identity in our <p 2 x 2 model; 



lim F 



1 

'2 



F(0) <p(x)<p(y) - F_(A) - {<p{x)<p{y) - (p(x)uj(y)) + F + (A) 

<p(x)<p(v) - \(vM*My) + <p(?)vMv)) + (m 2 + a x 2 M*M2/) 



A = ^/-V 2 + M 2 + \ x 2 , 

F(cj) + F(-cj) -2F(0) 



(CO 



F(u) - F(-cu) 



£0=2/0 
(2.40) 

(2.41) 
(2.42) 



2u 2 ' 2u 

Here we assumed, for simplicity and without losing practical utility for our prob- 
lem, that x 2 is independent of the spatial coordinate. For our purpose we may take 
(x 2 ) = T 2 /12 for this constant value, to give the temperature dependent cp mass 
M 2 + AT 2 /12 in thermal medium. We shall have more to say on the temperature 
dependent mass later. Furthermore, in the homogeneous thermal medium of our 
interest one may take the spatial Fourier component, 

J d 3 xe rk -*( V (x,x M0,y )} = (tf(k,x )<p(k,y )), (2.43) 
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thereby replacing V 2 by the relevant momentum, to get 



A = u k (T) = Jfc + M* + A T 2. 



(2.44) 



For a notational simplicity we often omit the temperature dependence of the (p mass, 
simply using the notation Uk for u k (T). 

Using the isotropy of medium, one has for the odd power terms, 

2 d d —> — » 

!im — (- — ) (<p ] (k , x )ip(k , y )) 

x ^ y + 2 dx dy 



-i , d 



d , 1 



lim 

xo ^ y + 2 dx dy 2 



4 v 



dfl(k,x ) 
dt 



<p(k,x ] 



(ip^(k,x )ip(k,y ) + 0(k,x o )0t(k,y o )) 
d0(k,x o ) 



+ 



dt 



ft(k,x ) 



) 



(2.45) 



using the equal time canonical commutator 

d(p{k , t) 



dt 



,0(-k',t) 



i5 3 (k-k>) 



(2.46) 



Even power terms, on the other hand, contain the modified occupation number, 



1 



dip^(k , t) dip(k , t) 



'<* • «> + 2 = 2^T) < "'-"a-' + ^(T) , , *) > ■ (2-47) 

Using this operator relation for 

2r x (A; + k' ,k + k') 



1 



-P(k +k' ) 



(2.48) 



one obtains the spectral weight for the Hartree model; the main result of this section, 
eqs. Q2.34j ) — Q2.35|) . Even and odd power terms have been lumped according to 



FijjJy) + F(-Uk') ( f(k , t) + ]\_ F(u} k ,)-F(-u k ,) 1 
2uj k > \ 2/ 2u k > 2 

F(w fc f(k' , t) + F(- w fc (1 + /(A;' , t)) ) , 



2cji./ 



(2.49) 



except the first power term F(0) giving the third term in the right side of eq.( [2.34j ). 
We shall later show that the combination of / and 1 + / correctly describes the 
destructive and the creative processes, including the stimulated emission effect for 
bosons. 
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It is important to realize that the spectral function in the reduced Hartree 
model, r(k ,t) of eq.( [2.34| ), cannot be written in terms of the distribution function 
f(k,t) = (aJ(k ,t)a{k ,t)) of (p particle alone, and contains another combination of 
the momentum and the coordinate, v(k) = {p\/{2uj^) — u k ql/2) , which vanishes 
for a system of a set of independent harmonic oscillators. To this extent the particle 
picture does not completely hold in a full quantum theory such as ours. In the next 
section we shall derive a self-consistent set of equations that determine the particle 
distribution function under the assumption of a slow variation of the particle number 
density. The v term becomes a higher order effect in A 2 . This way one can derive 
a useful quantum kinetic equation. It is convenient, with this in mind, to separate 
the contribution to the spectral function without the v term of ( 2.35Q , 



r J k t ) - 2 f d " k ' I r x(ko + ^k>,k + k>) g 

As a minor note, we point out that the relation f(k',t) = f(\k'\ ,t) following 
from the isotropy was assumed in the derivation above. We also ignored a possible 
tadpole type of contribution, which will be discussed shortly. 

The same operator equation, eq.fl2.40p, when applied to 



F(k ) = e iko(xo-yo) e -iHx-y) ? ( 2 _ 51 ) 

gives the full propagator in terms of the distribution function, 

-iG(x,y) = J {2 t k 2uk e^ 3 -® (e**<*>-*» f(k,t) 

+ e- iuJk{x °- yo) {1 + f{k,t)) -2v{k,t)) , (2.52) 



where u>k = y k 2 + M 2 + \\ 2 with a constant x 2 , an d t = (x + yo)/2. Again, the 
full propagator contains v in addition to the distribution function. 

In the previous discussion we ignored a possibility of a tadpole type of contraction. 
For instance, in the correlator of 

( V (IM2) f dxdy( V 2 (x)-^ 2 (x))(a(x-y) V 2 (y)- a *(x-y) V ' 2 (y))), (2.53) 
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one may contract both field operators at 1 and 2 with those of ip 2 at the same 
spacetime point x. This would lead to a new kernel of the form, oc 8(x — y) (ip 2 (0)) . 
More specifically, one allows a local term of the form, 

8(3 4 (x-y)=8(x-y)( V 2 (0)} f dya^x-y). (2.54) 

Jx >yo 

It is then easy to show that 

5(3 4 (x) = -i\ 2 5(x)4(ip 2 (0)) I" dko ^ 0,0 ^ coth^, (2.55) 

J 2m Ko 2 

with m the \ mass. There is another tadpole term to even lower order A from 



2 

which gives 



^M^(2) / dx V \x) X \x)), (2.56) 



5(3 2 (x) = 1 8{x) X ( X 2 (0)) = 1 8(x) X I - coth ^ . (2.57) 



(2tt)32^ 

These tadpoles, ( |2.55| ) and Q2.57| ), are shown in Fig.l. 

Both of these 8f3i(x) here are local (oc 8{x)) and are purely imaginary. These 
contain the mass and the coupling constant renormalization, the (f 4 coupling. Thus, 
the infinite part of these is eliminated by their respective counter terms. The remain- 
ing, temperature dependent finite terms give physical effects of finite temperature 
correction to the mass and the coupling constant. For discussion of order A 2 effects 
in the kinetic equation, only the temperature dependent mass contributes in the 
Hartree model, changing the renormalized mass to 

M 2 -> M 2 (T) = M 2 + ^T 2 . (2.58) 

The coefficient ^ was calculated from the vacuum subtracted A (x 2 (0)) — (x 2 (0))t=o 
above, and precisely coincides with the result of finite temperature field theory |15 |. 
The temperature dependent (p 4 coupling of order A 2 gives a higher order term to our 
subsequent kinetic equation, hence may effectively be ignored. 

The temperature dependent mass has the structure of the tadpole term, 

-i\{ X 2 )^\x)- V '\x)) (2.59) 

in the influence functional. In the usual operator formalism at zero temperature 
this type of term (x 2 ) is absent by the normal ordering, but in thermal environment 
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(x 2 ) has a temperature dependence, hence has a physical meaning. After renormal- 
ization one may write the temeperature dependent mass term as \M 2 (T) ip 2 in 
the Lagrangian density. With this new term, one should replace the previous a(x), 
eq.(£TUp by 

a{x) = \ 2 (trT[ X 2 (x)x 2 (0)] P/3 -(trxW) • (2-60) 

With this new contribution included, one can forget about the tadpole term provided 
that one uses the temperature dependent mass, M(T), instead of the renormalized 
mass parameter. We shall come back to further aspects of the renormalization in 
Appendix B after we introduce the slow variation approximation in the next section. 

The two-body spectral weight r x (k) for a(x) is calculable from the analytically 
continued imaginary-time thermal Green's function Hl2l , fl3fl . We first note the 



oddness, r x (—u,k) = — r x (uj,k), hence r x (— k) = —r x (k) combined with the 
isotropy of space. The spectral weight is given by a discontinuity along the real axis 
of the energy u = ko corresponding to two x particle states in thermal medium. Since 
the kinematics is modified in thermal medium from that in vacuum, the relevant 
expressions are different, depending on relation between k and \k\. For both k > 
k 2 + Am 2 and \k\ > k > §, 



A 2 / / Am 2 , 2 1 — p-Pw+ \ 



kn k / 4m 2 



where k = \k\. The formula is better understood |14| if one writes this separately in 
the respective kinematic regions, using the thermal distribution function /th(^) = 
l/(e^ - I)] 



for to > v k 2 + 4m 2 

A 2 r 



r x (" , k) = [ + dE ( ( 1 + f th (E) ) ( 1 + f th (u - E) ) - f th (E) f th (u - E) ) , 

(2.63) 



whereas for \k\ > k > 



2 



r x (co,k) = ^- J ^ dE(f th (u J )(l + f th (u J + E))-(l + f th (E))f th (u J + E)) . 

(2.64) 
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Thus, the indivisual discontinuity has a one-to-one correspondence to a physical 
process such as something «-> \X an d something + x ^ X- An example of the 
two-body spectral function r x is given in Fig.2. 

The spectral weight r (uj,p,t) in the Hartree model given by ( 2.5(Jp takes a 



particularly simple form if one uses the energy conservation for the the thermal 
factor, f^(u) + 1 = . For instance, near the mass shell, u ~ ±co> p , the relevant 
processes are the annihilation and the scattering process for uj « uj p and their inverse 
processes for uj = — uj p ; explicitly 

r (uj p ,p) = 2A 2 J dp' J dk\ J dk 2 

■ ( ( 1 + f th (h) ) ( 1 + f th (k 2 ) ) f( P ') 5(p + p'-k 1 - k 2 ) 

+ 2f th (h) ( 1 + f th (k 2 ) ) ( 1 + fip') )6(p + k 1 -p'-k 2 )), (2.65) 
r (—uj p ,p) = 2A 2 J dp' J dki J dk 2 

■ ( /th(Aa) U{k 2 ) ( 1 + f{p') ) 5(h + k 2 -p-p') 

+ 2f th (k 2 ) ( 1 + f th (h) ) f(p') 5(p' + k 2 -p-h)). (2.66) 

Four processes appearing in these equations are depicted in Fig.3. Here 

/th ^) = e ^l_ i ' °° k = ^ k2 + m2 ■ ( 2 - 67 ) 
We used an abbreviated notation for the phase space integral; 

*<■■■> = / (^r k -»' etc " (2 ' 68 > 

and S(p + p' — k\ — A; 2 ) = (27r) 3 5 4 (p + p' — k\ — A; 2 ) , etc. When multiplied by 
fip) and 1 + f{p), these give the on-shell destruction and the production rate of the 
momentum mode p with a suitable normalization. 



3 Slow variation and self-consistent equation 

Once the Hartree approximation is defined with the kernel f3(x ,y) in the influence 
functional given, one may proceed as follows. We have in mind a physical situation 
in which one can clearly separate two time scales; the microscopic time scale for the 
quantum behavior governed by the Hamiltonian, and a macroscopic time scale for 
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the change of the (p particle distribution function f(k,t). The separation naturally 
occurs for a small system within a dilute thermal medium at very low temperatures, 
where the system itself is not far from equilibrium. Under this circumstance one may 
consider a slowly varying f(k ,t) in the time range, t 3> microscopic relaxation time. 
We thus take, in considering the short time variation, the time t in the distribution 
function as a fixed constant. It is then straightforward to work out consequences of 
the short time dynamics, since the truncated Hartree model is formally equivalent 
to the harmonic oscillator model, known to be completely solvable ||, M. Since 
the relaxation rate towards the equilibrium depends on the distribution function at 
that time, this consideration leads to a self-consistent equation for the the quasi- 
stationary distribution function. 

Separation of the two time scales to distinguish the slow and the fast process 
also appears in other approaches fl|], || when the quantum Boltzmann equation is 
derived in different contexts. 

We shall briefly summarize the main point of how the exactly solvable model is 
used in our Hartree model. One first recalls that one can discuss each Fourier k mode 
separately. The Fourier mode kernel f3(k,t) in the influence functional is related to 
the kernel given at the relative coordinate x by 

P@,x Q ) = J ^- 3 (3(k } x )e rk -\ (3.1) 

For the sake of clarity we sometimes omit the mode index k unless confusion occurs. 
One notes that the path integral for the p system in the influence functional has an 
exponent of the form, 

l -^dr(i(r)X(r)-uj 2 (T)C(r)X(T)) 

- f tf dr [ T ds ( £(t) p R (r - s) £(s) + i £(r) (3 t (t - s) X(s) ) . (3.2) 

Here uj 2 (T) = k 2 + M 2 (T) includes the temperature dependent mass given by 
eq-(P38|)- 

The time Xq in eq. fl3.lp is a common CM time in the original f3(x ,y). Since there 
is only one common time, one has the time translation invariance for /3, which was 
absent in the original (3{x ,y). 

A remarkable feature of this exponent ( |3.2|) is that it is linear in the variable 
X = ip + ip', hence its path integration yields a trivial delta function. This gives an 
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equation for the semiclassical path for the £(= <p — <p') variable; 

+ co 2 (T) £(r) + 2 jf* rf S - r) = . (3.3) 

We set hereafter tj = and tf = t. This integro-differential equation can be solved 
by the standard technique of the Laplace transform, as described in ||. Its solution 
£d is given by 

UT)=£iV(T)+£fZ(T), (3.4) 

»M = ^, 2W ^( t - T ,-^ip, (3,) 

o(r) = / dp = — , (3.6) 

yV 7 27Ti ipo-ioo p 2 + cu 2 (T) + 2/3(p) 

/3(p)=/ dre-rrpfr), (3.7) 







with p > 0. Note that £(0) = 1. 
Applied to our specific model, 

= _ r iu ?ig4, (3.8) 

jo p^ + ClT 

r±(a> , fc) = r(a> , k) ± r(— a; , k) , (3.9) 
for the momentum mode k. The basic function g(r) for this mode is then 

9(k,r) = — / dze- izT F(z,k), (3.10) 

Z7T J-oo+i0+ 



z 2 — C<J 2 



(3.11) 



We deformed the contour of integration as p — > —iz'm the Laplace inverted formula. 
Both g(r) and its time derivative g{r) obey a related equation to (|3.3|) ; 

-j4 + ^ 2 ( t ) x ( t ) + 2 r ^ - s ) x ( s ) = ° > ( 3 - 12 ) 

with different boundary conditions, g(0) = 0, g(0) = 1 . 

Important quantities for the influence functional in the Hartree model are two 
real functions 



f3ji(k ,t) = / duor + (uj , k) cos(uit) , (3.13) 

JO 

_ /*oo 

/3Ak,t) = - dujr_(uj ,k) sm(ut) . (3.14) 

JO 
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Different combinations r± appear, because these integrals originally defined in — oo < 
oo < oo have definite parities; Pr(Pi) is even (odd) in t. These $ depend on /(A;) 
of the particle distribution via r±(u , k). 

For comparison we write the corresponding quantity (3^ in the unstable particle 
decay described by the Lagrangian density, A ipx 2 (A having a mass dimension in this 
case) ; 

(3ft\k,t) = ( doo coth — r x (oo , k) cos(otf) , (3.15) 
J o 2 

(3f\k,t) = - f°° dcor x {u,k) sin (cot) , (3.16) 



with r x given by eq. Q2.61l ). Thus, in this case r±^ that defines (3^ \ satisfies r+\oo , k) = 
coth^r { *\u,k). 

On the other hand, it turns out that in the present annihilation-scattering prob- 
lem a similar relation holds only on the mass shell and only when one takes the 
thermal distribution / t h for f(k); 

r + (oo k , k)th = coth r -(^fc > fc)th , (3.17) 

which further reduces to the detailed balance relation, 

r(uo k , k) th f th (k) = r(- uo k , k) th ( 1 + f th (k) ) . (3.18) 

A deep reason why this relation does not hold for the more general non-equilibrium 
case in our annihilation-scattering problem is that r±(u,k) are functionals of the 
non-equilibrium (p distribution and in general, f(k) ^ fth(k). In the simple case of 
the unstable particle decay into two thermal particles the corresponding spectrum 
is independent of the distribution function of the decaying particle. We further note 
that distinction of two types of the spectrum r± is crucial to obtain the correct form 
of the Boltzmann equation in the next section. 

Under a certain condition the function F(z ,k) is analytic except on the real 
axis where there is a branch cut singularity as shown in Fig. 4. As seen by taking 
the imaginary part of the defining equation ( |3.11| ), this analyticity holds under the 
condition, 

r_(oo ,k) = r(oo ,k) - r(—u,k) > 0. (3.19) 

This condition for the analytic property of F(z , k), when evaluated on the mass shell 
oo = oo k , physically means that the destructive process of <p given by the first term 
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r(u k , k) dominates over the production process given by the second term r(—u k ,k), 
a situation we are practically interested in. Furthermore, in the weak coupling 
limit the above condition needs to be obeyed only near uo = u k , since the off-shell 
contribution in the weak coupling limit is negligible in determining the analyticity. 
The analytic extention of F(z , k) has two simple poles in the second Riemann sheet, 
approximately at 

z = ± oj k (T) - i . (3.20) 

2u k 

The poles in the second sheet are close to the real axis in the weak coupling limit, 
A 2 — > 0. When pole terms dominate in the integral, the exponential decay law 
follows; 

^V),^!^, (3. 21 ) 
T k = - r ~ (uJk ' k) = - ( r (a;* , k) - r{-u> k , k) ) . (3.22) 

LO k U k V ' 

As is clear in ( [2.65 ), ( 2.66 ), the decay rate T k has both annihilation and scattering 
contribution along with their inverse processes. 

On the other hand, when r_(u k ,k)<0, these poles appear in the first Riemann 
sheet. The physical situation here is the dominance of the production process, which 
is not of our immediate concern. It however appears that this case can be dealt with 
analogously to the above case of r_ > 0. 

We now turn to derivation of the self-consistent equation for the distribution 
function. Recall first for each Fourier mode k, 

f{k,t) = (a}(k,t)a(k,t)) . (3.23) 

Here the creation (a') and the annihilation (a) operators of tp particle are Heisenberg 
operators such as e tHt a e~ lHt with H the total Hamiltonian including the system- 
environment interaction. In the harmonic oscillator basis which is an essential ele- 
ment of the plane wave decomposition of the field operator, this is equal to 

/(M) = (^j-Pk(t)p- k (t) +Y<lk(t)q- k (t) - ~> , (3.24) 

where the coordinate and the momentum, q k ,p k , are identified using the field de- 
composition into the plane wave; in the previous notation, 

q k (t) = 0(k ,t) , Pk {t) = (3.25) 
<p(k,t)= [ d 3 x V {x ,t)e~ in - s . (3.26) 
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The subtracted factor | in ( |3.24| ) is the well known contribution from vacuum fluc- 
tuation. For simplicity we subsequently use the notation, for example q\ to mean 
QkQ-k ■ 

We mention here an ambiguity for the choice of the reference energy uj^ to define 
the occupation number. We used here the "free" part of the oscillator energy Uk = 
\Jk 2 + M 2 . In the presence of the interaction with thermal environment the use 
of this unobservable energy is however dubious. Indeed, we confirmed, as will be 
explained in Appendix B, that this choice would lead to an unacceptable result of the 
absence of the equilibrium distribution at very low temperatures. As will be discussed 
later, the proper choice of the reference energy turns out to be the renormalized 
energy including 0[X 2 } correction of the finite self-energy shift in medium, eq( p.37p . 

We need to sharpen the concept of various time scales in computing the correlator. 
Consider the correlator at different times such as (q(t + t/2)q(t — t/2)) , and 
suppose that to ^ t, all the times measured from some specified initial time. We 
then let to 3> several x relaxation time, and vary the relative time t in the range 
from to several x the relaxation time. Under this circumstance we approximate 
(q(t + t/2)q(t — t/2)) « (q 2 (t )) . We thus need the coincident limit at large times 

relaxation time). 

Some means to compute the correlator at the coincident time such as (q 2 (t)) 
becomes necessary. This or a more general multi-time correlator such as 
(9(^1)9(^2) " " " q(t n )) can be computed with the machinery of the generating func- 
tional. Instead of being much involved in technical details, we shall give only a 
general idea here and relegate all technical points of the generating functional to 
Appendix A. An alternative method of computation is to use the exact operator 
solution for the harmonic system, as is done in |J. 

Calculation of the Green's function in the path integral approach is performed 
by introducing a coupling of external source terms of the form, 

j(r)q(r) + l(r)p(r)~f(r)q'(r)~l'(r)p'(r) = ± ( S£ + DjX + Sfo + D lPx ) , (3.27) 

where Sj = j + f , Dj = j — f , px ,§ = p ± p' , etc. Functional differentiation with 
respect to the sources, Si , -Dj then gives various combination of correlators, which in 
turn gives necessary two point correlators, (g(l)g(2)) , (p(l)p(2)) , (g(l)p(2)) . 
The result is 

(q 2 k (t))= / dur + (u;,k)\h(u,k,t)\ 2 
Jo 
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+ g 2 {k , t) p\ +g 2 (k, t) qf + g{k, t)g(k , t) Vi q { + q iPi , (3.28) 

POO _ 

(pt(t))= / dur+(u,k)\k(u,k,t)\ 2 
Jo 

+ g 2 (k,t)pj+g 2 (k ,t)~q? + g(k ,t) 9 {k ,t) p iqi + q iPi , (3.29) 

1 roo , _ _ . 

- (qk(t)Pk(t) +Pk(t)q k (t) ) = dur + (u,k) 3l[h(u ,k ,t)k*(u ,k ,t) J 

+ g{k ,t)g(k ,t)p! + g(k ,t) 9 (k,t)q[+ {g 2 (k ,t) + g(k ,t) g (k,t)) - mi + QiPi, 

(3.30) 

h{u,k,t) = ( drg(k,r)e- iujT , (3.31) 
Jo 

k(u,k,t)= f dTg(k 1 T)e- iuJT . (3.32) 



The coincident time limit in these formulas may be understood only with the con- 
dition, | relative time| <C t , thus the relative time can be as large as the relaxation 
time. Quantities such as qf are the values averaged over the ensemble at a specified 
initial time. We have chosen the initial ensemble such that ql = and pi = . Note 
that the relevant kernel for the correlator is /3ii(k,t), as seen from ( |A.19| ), hence 
the corresponding spectral combination is r + (u , k) appearing in the to integral here 
instead of r_ . 

We now coarse-grain the short time dynamics to derive the self-consistent equa- 
tion for the equilibrium occupation number. We first note that except at very early 
and very late times the simple exponential decay law for g(k ,t) is an excellent ap- 
proximation. This intermediate epoch of the exponential decay law is ideal to obtain 
the coarse-grained behavior. A sacrifice resulting from the replacement by the in- 
termediate exponential law is that the very early quantum behavior is lost. We lose 
nothing, however, in the late time behavior, because the power law behavior present 
at much later times is not realized due to a large multiplicity of actual reactions; 
indeed, the relaxation rate Tk which will be precisely defined later depends on the 
time weakly via the slowly varying distribution function f(k,t). 

We observe that in the formulas for the dynamical variables, eq.( |3.28| ) — ( |3.30| ), 
the initial value dependence disappears with the exponentially decaying g{k , t). We 
then use the limiting behavior of h(u , k , t) and k(u , k , t) as t —>■ oo, to actually mean 
t ^> relaxation time. From the expression of g(k,t), eq. (|3.10|) , and the definition of 
these functions in ( p.31 ), ( 3.32 ) we find that 

\h(u ,k,oo)\ 2 w — , (3.33) 

{u 2 -Cjl{T)Y + ^T_{u,k)f' 
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u 2 

\k(to,k, oo)\ 2 w — , (3.34) 

(^ 2 -^(T)) 2 + (7rr_(o;,A;)) 2 

c^(T) = £ 2 + M 2 + A r 2 + n(a; fe , fc) , (3.35) 
U{u,k) = -V I" dJ lAAAl . (3.36) 

A possible infinity in the proper self-energy H(uJk , k) is removed by a mass counter 
term such that the combination M 2 + Tl(u>k , k) is made finite by renormalization; 
after the wave function renormalization as explained in Appendix B, 

u>l(T) - (uf) 2 = k 2 + M 2 R + ^T 2 + 6H{u k , k) , (3.37) 



where Mr is the finite renormalized mass and 5U(uj k , k) is defined in eq.( p.ll 
The self-consistent equation for the stationary value is then 

f cq (k) + -= T d"(~ + — ) r+(h;,g) (3-38) 

J qK) 2 Jo K 2 2u k } (uj 2 -Cj 2 {T)) 2 + {ixr_{uj,k)) 2 1 ' 

v eq (k) = r *< ■ (3-39) 

q Jo 2uo k (u 2 -u 2 {T)) 2 + {7cr4uJ,k)) 2 



r±(u,k)=r(w,k)±r(-u,k), (3.40) 

r(u,k) = 2 I , f*' [ r ^ + ^f + k') g 
1 (2n) 3 2u k > I l- e -^+«*') ,/eqV 



+ 



In this self-consistency equation the two-body kernel r x is a given function. The 
function v eq (k) is the stationary value of v(k , t) given by 

<£(*)>. (3-42) 
and may be considered as a functional of f C q{k) when one first solves eq. ( |3.39| ) for 

Veq- 



The set of self-consistent equations, eqs.( p.38| ) — fl3.41|) , along with the definition 



of r x (|2.61| ), does not depend on the initial (p state. We assume that there exists 
a unique solution to the self-consistent equation, when a proper renormalization is 
performed. We also anticipate and later deomonstrate that the high temperature 
limit of this equation gives the familiar distribution function l/(e^ fc — 1). 
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A simplified computation is possible in the weak coupling limit. In the limit 
of A — > 0, one can ignore the u dependence of r_(cj ,k) in the denominator of the 
Breit-Wigner function since it is narrowly peaked around u = u k . By separating the 
narrow width contribution, one has 

Uk) = r{ ~ Uk ^ ] +5~Uk), (3.43) 

of cq {k) = - — / dw ~ ( rrw i r2/ o 3 - 44 

Auj k J-oo [u - u) k {T)Y + Iy4 

with Y k = 7r r_ {uj k , fc) / ui k . We have used r + (—u,k) = r + (u ,k). The tilded 5f eq 
contains terms to be renormalized away by subtraction, in contrast to the finite Sf eq 
later defined after renormalization. Note that both 5f cq and Dqq £1X6 of the same 
coupling order, 0[A 2 ], since the on-shell term of O[X ] is absent. It is then easy to 
see that the quantity v eq gives an even higher order 0[A 4 ] correction to r±, and one 
can forget about v eq altogether to 0[A 2 ]. 

The self-consistent equation in the present form is not particularly illuminat- 
ing, since cancellation occurs among scattering terms in the on-shell part when the 
distribution function is integrated over momenta, namely in the quantity, 



/ 



" Sk Tt ( U k) - r -^4 I • (3.45) 



In the next section we shall directly work out the stationary number density inte- 
grated over momenta. In the rest of this section we shall focus on the off-shell part 

Sfeq- 

To further simplify the off-shell contribution, we use a sum rule that results from 
a consistency condition; from the equality, g(0) = 1 , 



JO 



ujr_(u,k) .„ 
du — = 1 . 3.46 

(cui-Cu k (T)iy + (nr_(tu,k)y 



Using 



-oo iuj - + m ~ ^ ' (3 - 47) 

— # — # 

u)r(u,k) 1 f°° , r(u,k) 

' du ■ 



f°° , ur_(uj,k) 1 f a 

/ du =; — pa - /rTnX / 



(u 2 - u k {T) 2 ) 2 + (7rr_K , k)) 2 ' " 4u; fc (T) J-oo " (u - u k {T)) 2 + T 2 /4 ' 

(3.48) 
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with T k = -Kr_(oj k , k)/uj k {T) w nr_(uj k ,k)/u k , the consistency integral ( |3.46| ) may 
be rewritten as 

- r_(o;,fc) -r_(^ fc ,fc) 

oo^(.-^(T)) 2 + rl/4-°- (3 ' 49) 

We have neglected the a> dependence of r_ (u , fc) in the denominator in deriving this 
equation. 

Since r±(u) = r(ui) ± r(— cu) , the off-shell equilibrium distribution becomes 

8feq[k) = ~ / dw- , 3.50 

2w A 7-oo (a; - uj k (T)) 2 + Iy4 
K-«,fc =2 / ttt^tt- * ' , } f(k',t) 



(27r) 3 2a; fc / I e^-^'l - f 

+ rxil %Z k ' k !\ T i + ^ ( 1 + /(# ' ) + g(a; y - r x (|a; - u„\ , k - 9) f(k> , t) 

+ e(-uj-u k/ )r x (\uj + u k ,\,k + k')(l + f(k',t))) . (3.51) 

In the spectral function r(— , fc) here, one can disregard the term involving t> eq in 
eq. (|3.41|) , since it is of higher order, 0[A 4 ]. 

A region of u integral in the formula ( |3.44j ) for 5/ eq does not contribute; this is 
the region of \u — u k \ < u c , where uj c is the energy scale for which r(u , k) appreciably 
varies and in the weak coupling limit u c 3> T k . This means that the resulting 5f eq (k) 
is insensitive to the actual value of IV We may use this fact to replace T k by a small 
quantity T k which is taken independent of the distribution function f(k); 



LP ,k) - r(-u k , k) 
2u k J-oo (u-u k (T)) 2 + t 2 k /A 



suk) - ^ / r izZT £r; ■ (3-52) 



We shall specify T k shortly in ( |3.53| ) by a quantity even independent of the mo- 
mentum k. A great virture of this formula is that the off-shell part Sf eq of the 
self-consistency formula becomes a linear functional equation of f e(l (k), and the per- 
turbative treatment becomes transparant. This form of the equilibrium distribution 
is most convenient and is later used frequently. 

The universal width factor T k independent of / may be defined by using the 
value of the width T k calculated at / = 0. Since T k = nr^(uj k ,k)/u k is the on-shell 
value, it is dominated by the scattering contribution, T k = 7ir s (uj k ,k)/u k . Thus, 
it numerically has the thermal number density factor n t h ~ T 3 times the scattering 
cross section of order X 2 /M 2 . More explicitly, r*. becomes independent of the mode 
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k, hence is denoted by T with 



f A2T3 

T = ^AP- (3 - 53) 



In Appendix B we identify the renormalization term in the proper self-energy 
H(u , k). This gives constant counter terms to be added to the distribution function; 
the renormalized finite occupation number is 

fUk) = f cq (k) - ■ (3-54) 

The infinite counter term A , B cancells the corresponding infinity of the 0[X 2 } cor- 
rection in the distribution function. 

In Appendix D we give a detailed account of our computation of 

r-{uj k , k) 

To the leading order of T/M the result of this computation is summarized as a sum 
of a / dependent and a / independent term, 

5f eq (k)=f f (k) + f°(k). (3.56) 

As will be explained in Appendix D, the /—independent term has a dominant con- 
tribution from the inverse process XX ~^ PP, an d is computed as 

1 r d 3 k' 



f ° 2{k) ~ 2co k loo dU (u - u k y + 11/4 
r x (\u + ujy I , k + kl 

e /3\u}+u) k ,\ _ i 

2 r d 3 k' 1 



u k J (27r) 3 2cj fe / (u k + uj k >) 2 Jo 
C(2)A 2 T 2 f°° q 




(3.57) 



16tt 4 ku k Jo d<1 2q + ((2)T e<?/ 2T - 1 [u k + u k - q u k + u k+q ) ' (3 ' 58) 

The /—dependent term // is given in Appendix D, but actually is not necessary 
to write the time evolution equation for the number density. The use of the low 
temperature formula for r± for computation of 5f cq even at high temperatures is 
justified, since at higher temperatures the on-shell contribution dominates and the 

— * 

term 8f eq becomes irrelevant in f eq (k) which is dominated by the familiar / t h- 



28 



An example of the us— integrand ( [3.571 ) is shown in Fig. 5, where the total, the 
Planck, and the rest of the integrand are separately given. 

Physical processes that contribute to the important piece /| are predominantly 
inverse annihilation XX ~^ Wi an d 1 to 3 process, x ~^ X^P^P-, which gives a smaller 
fraction of the total number density. On the other hand, the term ff arises from 
the inverse scattering process. We shall show in the next section that the scattering- 
related term fj(k) is subdominant compared to f® (k) in determining the equilibrium 
number density. 

Although it is technically complicated to calculate the equilibrium distribution 
function / eq , one may regard this calculation as a self-consistent approximation for 
the quantity, 

/e q (fc) tr e _ pHtot , (3-59) 

where H tot is the total Hamiltonian including interaction between the (p system 
and the x environment. We explicitly demonstrated this relation for the exactly 
solvable harmonic oscillator model in ||, whose result is used for our self-consistent 
solution. Thus, the Gibbs formula is valid, while the ideal gas form of the distribution 
l/(e^ Wfc — 1) is changed. 



4 Quantum kinetic equation 



We go back to the time dependent expectation values, given in eq. (|3.28|) — ( |3.30| ) 
It is sometimes convenient to write a formula for the time derivative; 

^M = -P(M) 
at 

■(f(k,t)- [°° dwr + {uj,k){^\h{u;,k,i)\ 2 + ^-\k{u;,k,t)\ 2 
V Jo 2 2,00k 

1 

2u k 



+ r^^ r+ (u;,fc)(^|^,fc,t)| 2 + -^|^^,t)| 2 ), (4.1) 
Jo at 2 2uJk 



r(k,t) = 



• 2 



t: In (— 9 + — Pi + (it 9 + 7^— it + o K 9 + ~ ) QiPi + PiQi 
at \ 2 2oo k 2 2uJk 2 uo k 



(4.2) 
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Regarded as a time evolution equation, this equation has a remarkable property; 
the initial memory effect is confined to, and isolated by, the rate T(k,t). All other 
quantities in this equation are written in terms of those at the same local time, 
including the distribution function f(k,t) itself. 

We take a Markovian limit of this equation. The idea is as follows. The equation 
above, although exact in the slow variation limit of the Hartree model, has the 
initial memory effect. To retain the memory effect to late times may not make much 
sense when one wants to discuss an average behavior of the time evolution for a 
collective body of particles in a complex environment. After all, it is impossible to 
follow the time evolution for all particles in the ensemble. It may be possible to 
specify the initial data such as the ensemble-averaged qf for all momentum modes at 
some specific time, and the data may have a simple regulated form if one assumes a 
simple initial distribution function for (p. But, remember that in the problem of our 
interest there are processes occuring rapidly such as the scattering, distinct from the 
other simultaneously occuring slow process such as the annihilation, which is of our 
main interest. Under this circumstance different particles undergo different history 
of evolution, and after a while it is almost impossible to keep track of the updated 
initial data. 

One really does not care much about the details of the complicated time history. 
We only care about a global and slow change for a few key quantities. Fortunately, 
in many physical situations one has a clear separation of at least two time scales; 
one due to an elementary quantum behavior and the other for the bulk behavior. 
With this in mind it is much more sensible to erase the initial memory effect and to 
coarse-grain the time evolution averaging out the fast oscillatory behavior. A nice 
feature of the Hartree model is that if one eliminates the initial time dependence 
in T(k,t), a simple Markovian description of the transport phenomenon becomes 
possible. 

In the weak coupling limit the exponential law g(k,t) ~ sin(a;& t) e~ Fkt ^ 2 /cu k 
gives a constant decay rate, 



T k = — r_(u k ,k) = — (r(u k ,k) -r(-u k ,k)) , (4.3) 
u k u k v ' 

after the time average. Except at very early times the initial memory effect is almost 
completely erased. Via the spectral function this T k depends on the (p distribution 
function. As seen from eqs. (|2.65|) , (|2.66 ), this combination r_ is the rate of decrease 



30 



of heavy particles due to the annihilation, the scattering and their inverse processes. 
Our Markovian evolution equation is then 

#££ = -r u (f(k,t)-Uk)), (4.4) 

Uk) = r{ ~ U ""$ l + SUk). (4.5) 

Here the off-shell contribution 5f eq {k) is given by ( |3.56|) — ( |3.58| ). Although not 
written explicitly, the distribution function used to compute T k and f eq (k) should be 
the instantaneous one, f(k,t), with the same common time t as in the left side. The 
Markovian evolution equation for the momentum distribution function thus derived 
is the most fundamental result in the present work. 

We note that the same Markovian approximation, when applied to the model of 
unstable particle decay, gives the kinetic equation identical to eq.(fOJ), except that 
Tk in that case is the decay rate on the mass shell and that r±{u , k) in f eq {k) should 
be replaced by the appropriate spectral weight for the decay process such as r x (u , k); 
more precisely r + — > r x coth 4^ , r_ — > r x . We give in Appendix C the expression 
and its actual value of the stationary distribution f eq (k) for the boson decay model. 

We shall recapitulate the problem related to v(k). The Markovian equation 
derived above does contain the spectral combination r±, which is written using the 
distribution function f(k,t) and and also v(k,t). This new function describes a 
deviation from the simple particle picture in the field theory; in terms of the harmonic 
coordinate and its conjugate momentum, 

v(k,t) = (^-pt(t)-^ q l(t)}. (4.6) 

This combination of dynamical variables obeys another time evolution equation dif- 
ferent from the distribution function. The coarse-graining, under the same slow 
variation approximation applied previously, gives 

v$,t)= r du^i r -±^ (4.7) 



The initial memory term for the combination ( |4.q ) has the fast oscillatory behavior, 
which vanishes on the average. This behavior is different from that of the occupa- 
tion number which has a slowly decaying component, non- vanishing after the time 



average. In this sense it is best to use the expression fl4.7|) for v(k , t) and the defin- 
ing equation of r± containing both f(k,t) and v(k,t), as a functional equation to 



31 



determine v(k,t) in terms of f(k,t). This way we arrive at a closed form of time 
evolution equation for the distribution function. In the Markovian time evolution 
equation the quantity v(k,t) appears via the spectral function r(±u ,k), which we 
now show to be of negligible higher order. 

In the weak coupling limit of our main concern, both r± has an overall small 
coupling, and the integral ( |4.7| ), that excludes contribution from the Breit-Wigner 
region of to ~ Uk by the factor u 2 — uj\, gives v = 0[X 2 }, hence gives 0[A 4 ] contribution 
to the spectral function r, as seen from ( |2.34| ). We shall thus take v vanishing from 
now on. 

With the vanishing v, the usual Boltzmann equation follows when one approx- 
imates the energy integral for / eq by the pole term, namely r(—Uk , k)/r^(ujk , k) . 
This narrow width approximation gives the evolution equation, 

= -r k ( f(k,t) „ ) . {u 

at \ r(uk , k) - r(-u k , k) 



A 2 j dk' J dh J dk 2 



One may use the formula ( f4.3|) for and write this as 

VihA = -JL f r (u; k ,k)f(k,t)-r(-u; k ,k)(l + f(k,t))) . (4.9) 
dt ujk 

Again, the instantaneous ip distribution function is taken in computing r(±Uk,k) 
of the right hand side. In view of the on-shell relation (|2.65 ), ( 2.66|) for r(±Uk , k), 
this is equivalent to the Boltzmann equation for the present annihilation-scattering 
problem in thermal medium; 

df(k,t) 
dt 

( ( 1 + /th(A;i) ) ( 1 + U(k 2 ) ) f(k') f(k) 5(k + k'-h- k 2 ) 
+ 2/ th (^) ( 1 + f th (h) ) ( 1 + f(k') ) f(k) 5(k + kx — k' — k 2 ) 

- U(ki) U(k 2 ) ( 1 + f(k') )( 1 + f(k) ) S(k + k'-k 1 - k 2 ) 

- 2f th (k 2 ) ( 1 + fth(ki) ) /(AO ( 1 + f{k) ) S(k + k 1 -k'-k 2 )) , (4.10) 

where all relevant processes are explicitly written. We thus find that our Markovian 
Hartree model gives a foundation to derivation of the Boltzmann equation under the 
narrow resonance approximation. 

Note that the thermally averaged rate in the right side of the Boltzmann equation 
( |4.10| ) has separate contributions of the annihilation and its inverse process charac- 
terized by 5(k + k' — k\ — k 2 ), and the scattering plus its inverse by S(k + k\ — k' — k 2 ). 
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There is no contribution from 1 to 3, and 3 to 1 processes such as (p <-> pxx m the 
Boltzmann equation. When one integrates over the particle momentum k to discuss 
the time evolution of the number density, the entire scattering contribution drops 
out. This is reasonable, because the scattering process does conserve the particle 
number, hence the scattering process does not cause the change of the (p particle 
number. However, the momentum distribution function changes by the scattering. 
Moreover, it is not immediately clear how the off-shell scattering term contributes 
within our crossing symmetric approach. We thus keep the scattering term for the 
time being, as it stands. 

The equilibrium distribution function is defined by setting df/dt = 0. We as- 
sume as before that this has the unique solution, which in the case of the Boltz- 
mann equation must agree with that of the thermal and the chemical equilibrium, 
f(k) = \/{er Uk — 1) . The concept of partial equilibrium may however be of some 
use. For instance, if the scattering process takes place much more frequently than 
the annihilation process, it may be useful to suppose that the energy exchange is 
equilibrated, but the species number change is not fast enough. In this case the 
thermal equilibrium, and not the chemical equilibrium, is reached with 

where fi is the chemical potential. Under this thermal equilibrium of a finite chemical 
potential all scattering- related terms in eq. ( }4. 10|) cancels each other since 



[r s (uj k , k)} , tf h (k) - [r s (- u k , k)} , ( 1 + /&(*) ) = , (4.12) 
and only the annihilation-related terms are to be retained. Here 

r s (uJk , k) = 2 A 2 J dk' J dki J dk2 

■U(h) ( 1 + f th (k 2 ) ) ( 1 + f(k') ) S(k + h-k'- k 2 ) , (4.13) 
r s (-cu k , k) = 4A 2 J dk' J dk x J dk 2 

■U(k 2 )(l + / th (A;i) ) f(k') S(k + k\ — k' — k 2 ) . (4.14) 

In this case the Boltzmann equation effectively describes the time evolution for the 
chemical potential. Note however that the thermal equilibrium distribution of a 
finite chemical potential is realized, assuming that the on-shell Boltzmann equation 
is a valid description of our problem without the off-shell contribution. 
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As will be made clear shortly, it is not easy to obtain a readily calculable form 
of the equilibrium distribution function at low temperatures in the annihilation- 
scattering problem, and we shall directly work out the integrated number density. 
On the other hand, for the unstable particle decay we have an approximate analytic 
result, eqs. flC.lO ) and (|C.13| ), for the distribution function, which is shown in Fig. 13. 
At low T < 0.1 M deviation from the Planck form becomes large, but the distribution 
is not described by the form ( j4.11| ) with a finite chemical potential. 

More generally, it is convenient to separate the on-shell term and write the rest 
of contribution, as is done in the preceeding section; 

Uk) = - r( -^; fc) + 5f eq (k) , (4.15) 
r(u k , k) - r(-u k , k) 

5f cq (k) = f f (k) + f° 2 (k), (4.16) 

with f = ■ The function f 2 is given in eq.( |3.58|) . Physically, the off-shell 

contribution 5f eq (k) in this formula consists of two terms; the first / dependent one 
// due to the inverse scattering process, and the second / independent one f 2 due 
to the inverse annihilation process, XX ~^ W > along with a small contribution from 
1 to 3 process, \ XW • 

For discussion of the time evolution, we shall be content with the integrated 
number density, 

n(t) = J ||,/(M), (4-17) 

and its evolution. In the discussion of the relic abundance of WIMP this integrated 
quantity is of prime interest in cosmology. As already noted, the evolution equation 
for the number density is simplified considering cancellation of the scattering terms 
in the on-shell part, 

dn r d 3 k 

~dt ~ " ~ J (2tt) 

where r^ nn (r 1 ^) is the rate keeping the annihilation (inverse annihilation and inverse 
scattering) term. The approximate form of the evolution equation for the number 
density at low temperatures (T M) is then given by 
dn r d 3 k n 



/ 7S5 ( r r /(m) -n™ uk,t)) , (4.i8) 



* iw J'.K,t)/W-iTf«t)), (4.19) 

r a (uj k , k) = 2 A 2 J dk' J dk x J dk 2 

•( 1 + f th (h) ) ( 1 + f th (k 2 ) ) f(k') S(k + k'-h- k 2 ) . (4.20) 
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The equilibrium number density is obtained by setting dn/dt = 0, hence is given 

by 

RHS = - J (r a (cu k ,k)f(k)-rr5Uk))=0. (4.21) 

A more explicit form of this equation is 2A 2 times 



dk I dk' J dh J dk 2 [(l + f th (h) ) ( 1 + f th (k 2 ) ) f(k') f(k) 5(k + k'-k 1 - k 2 ) 
( 1 + U(k 2 ) ) Sf eq (k) ( 1 + f(k') ) 5(k + h - k' - k 2 ) ] = , (4.22) 



where we used the abbreviated notation for the phase space integral ( |2.68| ). Thus, 
r™ v is given by the second term of eq.( f4.22[ ). 



Some sort of averaged cross sections are here in this equation when A 2 is mul- 
tiplied, and it would be useful to crudely estimate these. Let us assume that ip 
particles are non-relativistic and ignore the effect of stimulated emission. One has in 
the first term an averaged annihilation cross section a^v times the ip number density 
squared n 2 , where v is the relative velocity between two annihilating particles, and 
<j a v is the invariant cross section. Since the annihilation cross section slowly changes 
with energy at low temperatures, one may take the zero energy limit of the cross 
section, 

A 2 

16ttM 2 k j 

In the second term there are two contributions corresponding to the two / in 
Sf eq , ff and f 2 ■ One of these // is multiplied by the scattering cross section, 

A 2 

a s v s w - — — , (4.24) 

while the second f 2 is multiplied by o a v a . The first one is given by a s v s times the 
quantity of order n til ■ n[ , with n th = ^r-T 3 , hence, using ( [D. 10| ) , 

"•^li^ 2 ^' < 4 - 25 > 

The second one is given by a a v a times 

o cC(3)A 2 T 7 

n t h n° 2 ~ — — ~~z , 4.26 

2 192tt 5 M K ' 

with c ~ 0.27, eq.( |C.15| ). Thus, once the ip number density becomes 0[T 3 ], the 
second contribution is dominant by 0[MT 2 /n v ]. Equating this cr a f a n t h^2 to the 
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annihilation rate a a v a n 2 gives 



1 /cC(3) , T 



7T Z 



1927T 



A \ / — — T 

M 



0.0023 x A J-^T 3 
V M 



(4.27) 



We thus derived the equilibrium number density roughly of order 

~T 



10" d x A 




(4.28) 



at low temperatures which may become much larger than of order (MT) 3 / 2 e~ M l T 
determined from the Maxwell-Boltzmann distribution of zero chemical potential. 

This argument shows that with n° = 0[X 2 T 4 /M] the /—dependent off-shell 
contributions // of order n{ = 0[X 2 (T / M) a n^] are subdominant unless 



(4.29) 



Even the possibly largest case obtained numerically, eq.( p.!4l) , gives a ~ 1.35, thus 
confirming that the dominant off-shell contribution is n®, eq.(P-21|). 



5 Application to cosmology: relic abundance 

The cosmic expansion has a drastic effect on the annihilation process of heavy 
stable particles. The temperature of cosmic environment particles such as the \ 
particle in our toy model decreases with the scale factor, T w l/a(t) , which in 
turn results in less frequent reaction. This gives rise to a phenomenon called the 
freeze-out or the decoupling of the process. 

The freeze-out is described introducing a term of the cosmic expansion in the 
evolution equation. For the evolution of the distribution function f(k,t), the time 
derivative operator is modified to 

dt a dk" [ } 

Here a /a is the Hubble parameter H. When one integrates over the phase space to 
get the number density, the left hand side of the evolution equation becomes 

^ + 3Hn = -- - , (5.2) 
dt 
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assuming that the distribution function is sufficiently damped in the high momentum 
limit. The evolution equation is further simplifed by introducing the relative yield 
Y, 

Y=^r 3 , (5.3) 

since 

dY d n m o . dn . , 

This holds owing to the temperature-scale relation, T oc 1/a , hence H = —T /T . 
The approximate Markovian kinetic equation in the expanding universe is then 

dY 2A 2 r r r r 

— = - 2tt — J dk J dk' J dk x J dk 2 

( ( 1 + / th (A;i) ) ( 1 + f th (k 2 ) ) ( f(k') f(k) - f MB (k>) f MB (k) ) S(k + k'-h- k 2 ) 

- /th(fci) ( 1 + U(k 2 ) ) ( 1 + f(k') ) 5f cq (k) 5(k + h-k'-k 2 )), (5.5) 

where we introduced the Maxwell-Boltzmann distribution for the zero chemical po- 
tential, 

/mb (*) = e~ M l T exp(- ^) . (5.6) 

A very crude estimate of the freeze-out temperature goes as follows. One equates 
the equilibrium annihilation rate ~ a a v a ■ n v to the Hubble rate, 



T 2 /47T 3 N i— 

H = d , d=\— xil.ffiVN, 5.7) 

m p i V 45 

where N is the number of particle species contributing to the cosmic energy density. 
This argument, when applied for n v ^> umb, gives the freeze-out temperature, 

Tt » 700 x iV 1 / 3 (— f 3 ^ » 1.3 x 10 5 (^ ) 2 iV 1 / 3 (— ) 2 / 3 M 
m p i A 2 A m p i 

We used the unit of A, anticipating a strength of order the electromagnetic interaction 
A ps e 2 = 47tq;. The use of the off-shell formula for n v is usually justified for A = 0[e 2 ], 
since at low temperatures n v > n MB . 

The freeze-out yield is defined by Yf = (n^/T 3 ) T=T/ . We find for the range of 
parameters, A > 9.3 x 10~ 5 (M/GeV) - 32 , 10~ 3 GeV < M < 1 TeV 

Y f ps 0.06 x N 1/6 (— ) 1/3 ps 2.4 x 10~ 8 A^ 1/6 (M/Ge\/) 1/3 . (5.9) 
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Remarkably, this quantity is insensitive to the coupling constant A. Prior to the 
freeze-out epoch, this value Y does vary, but only gradually, since oc 
The freeze-out yield Yf is almost invariant in the rest of cosmic expansion, as will 
also be discussed in analytic estimate below. The present relic mass density is then 
estimated from 

p = MY f T*, (5.10) 
with T the present microwave temperature of ~ 3 K . Numerically, 

po « 4.1 x 10 4 iV 1 / 6 (^) 4/3 eVcm- 3 . (5.11) 

The closure mass density of order 2 x 10" 29 g cmT 3 pa 10 4 eV cm~ 3 requires that 
M < 0[lGeV], assuming N = 43/4. Thus, WIMP in the mass range far above 
1 GeV is excluded for the S-wave boson-pair annihilation model. 

Extension to the annihilation in a higher angular momentum state is of great in- 
terest, since LSP in SUSY models pair- annihilates with a large P-wave contribution 
||. This P-wave annihilation for LSP is related to the Majorana nature of LSP. It 
is beyond the scope of the present work to accurately calculate the annihilation rate 
in supersymmetric models, taking into account the Majorana nature and all con- 
tributing Feynman diagrams. It is however not too difficult to qualitatively estimate 
the P-wave annihilation rate, by simply taking into account the momentum, hence 
temperature dependence of various rates, (k 2 ) oc T. 

A more detailed behavior of the (p number density may be worked out by exam- 
ining 

(171 T 1 

- + Hn = - (a a v a ) (n 2 -5 (-)*>+ 1 T 6 - n 2 MB ) , (5.12) 

^ = 16^^' ' = - 27x W' (5 ' 13) 
n MB = (^) 3/2 e- M/T . (5.14) 

Our toy model gives the S-wave annihilation with p = 0, while p = 1 for the P-wave 
annihilation. We have extended the S-wave annihilation to the case of higher angular 
momentum without changing the effective coupling constant 5 relevant to our toy 
model. Equivalently, using the inverse temperature and the yield, 

" (Y 2 ~Y 2 ), (5.15) 



dx xp+ 2 y eq 
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X = T- v = ^dt- (5 - 17) 

The parameter rj is roughly the (on-shell) annihilation rate a a v a T 3 divided by the 
Hubble rate at the temperature equal to the particle mass, T = M. 

We plot in Fig. 6 — Fig. 8 a typical solution to the time evolution equation, ( |5.15|) ; 
Fig. 6 and Fig. 7 for the S-wave annihilation and Fig. 8 for the P-wave annihilation. 

The analytic estimate of the freeze-out temperature and the freeze-out yield JTTJ , 
which well reproduces the numerical estimate above, is as follows. One may consider 
with a good precision that the yield follows the equilibrium abundance Y e q until the 
freeze-out temperature. This temperature Tf is given by 

~dx^ = ~xY 2 ' q ' 
since after this epoch the inverse process is frozon and the yield follows 

dY_ _q_ 

dx x p+2 

Integration of this equation gives the final yield 



Y l . (5.19) 



p+i 

which agrees with 



Y(x) = ^ r-T- , (5.20) 

1 ' l-Vy-Mx-P- 1 -^ 



as T — > 0. Usually Yy is very small along with xj 1 1, and in this case Y ~ Yy 
after the freeze-out. 

In Fig. 9 — Fig. 11 we show the present WIMP mass density in the parameter 
space (M,A). The off-shell dominance region is shown in Fig. 9, along with the 
closure mass density, while contours of smaller mass densities are shown in Fig. 10. 
In Fig. 11 the P-wave case is shown. The excluded region due to the overclosure is 
larger for the P-wave annihilation than for the S-wave, with the same set of (M , A). 

For a given temperature T, there exists an upper bound on the heavy particle 
mass M max that can be produced in the equilibrium abundance n eq . In the Boltzmann 
equation approach this bound is roughly of order T, but our new contribution <5n eq 
substantially changes this value. There are two considerations to be taken here. The 
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first one is the condition on energetics; produced energy < thermal environment 
energy. From Mn cq < cT 4 with c a constant of order unity, 

M < M max , M max « 2 x 10 5 c 2 J . (5.22) 

This bound is not stringent, since M max /T 10 5 /A 2 can be quite large. 

The second, a more important constraint comes from the relaxation time. For 
the inverse process XX ~^ W to occur frequently, its rate r inv must be larger than 
the Hubble rate H. With 

r inv ~ (<r a v) n cq » 3 x 10~ 4 A 3 (-^) 7 / 2 M , (5.23) 

this consideration gives 

M max » 3 x 1(T 2 A 6/5 iV- 1 / 5 (^?) 2/5 T . (5.24) 

As T -> 0, M max -> like T 3 / 5 , but M max /T can become very large. 
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Appendix A Generating functional in the influ- 
ence functional method 



We explain the technique of the generating functional applied to the influence 
functional, taking the example of the exactly solvable model |§. When one wants to 
apply this method to the Hartree model of our annihilation-scattering problem, one 
should use relevant spectrum r± instead of r{u) below. 

We first introduce the source terms coupled to both the harmonic coordinate and 
the conjugate momentum as 

j(r) q (T) + l(r)p(r)-j'(r)q'(r)-l'(r)p'(r) = ^(Sjt + DjX + S& + D lPx ) , (A.l) 

where a convenient combination for the influence functional method is introduced 



Sj = J + J , Dj = j — f , Si = 1 + 1', Di = l-l' , p x ,f=p± p 



(A.2) 



In subsequent formulas we often omit the momentum index k for the sake of sim- 
plicity and discuss each Fourier mode separately. Funtional differentiation of the 
resulting density matrix pC?>0 with respect to these sources, when evaluated at the 
vanishing source, gives various combinations of correlators; for instance, 

5 2 



Sj{Tl)5j{T 2 ) 

5 2 

+ 



tr p 



(j ,1) 



(A.3) 



J j=0 ,1=0 



+ 



5 7 



SSfrdSSfa) 8S s {n)8Dfa) SDfrtfSfa) 
8 2 



+ 



(p{ti)p{t 2 )) 



tr p 



6 2 



J 3=0 ,1=0 

U ,0 



5l{n)5l(T 2 [ 



tr p 



(A.4) 
(A.5) 



J j,=0,Z=0 



Computation of the new density matrix p^ > l > under the action of the source is 
similar to the case without the source terms, because introduction of the source does 
not change the Gaussian nature of the Hartree model. The semiclassical £ equation 
and the the effective action is thus given by extending the analysis sketched in the 
text; 

Ur) = -£ f 9(t-T)+£ f g(t-T)+ f ds g{s - r) {D 3 - D^s) , (A.6) 
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U ,0 



iS. 



dr ds^ c x(r)/3 R (r - s)£ cl (s) 
Jo 

f dr (S s {t)Ut) + St^Ur) ) + % - (X f i f - 



2 Jo 

Here the change of the variable, — > is computed using 



with 



6 = "^W+^W+ / dsg(s)(D j - Di)(s 



1 Z" 00 

<#) = — / da;e-^F(a; + iO + ), 



dujr(uj) e 



-iuit 



F( I )-' = I »-^CT)-2 fi^M 



— C<J 2 



(A.7) 
(A.8) 

(A.9) 

(A.10) 
(A.ll) 

(A.12) 



For calculation of the correlator one convolutes the J^ '^ function ( |A.7|) above 
with an initial density matrix of tp system piiXi , £;) and traces out the final Xf and 
£j variables. Thus, it is convenient to take the trace with regard to the final variable, 
to get 



tr J™-. 
<M & - / drg(Dj - A) ) exp 



rfr / ds^ (r)/3 R (r - s)£ (s) 



+ 2~ i/ r (^(^o(r) + ^(r)e (r)) + ^ dr g{r){D 3 {r) - D^r)) 



dsiDj-Dt^gis-r) 



, (A.13) 
(A.14) 



As an example of the denisty matrix with the external source attached, we may 
work out the case for initial thermal state of temperature To = l//?o, to obtain 

-t 



trp 



exp 



f dr [ ds£ (r)f3 R (r - s)£ (s) 
o Jo 



+ -J o dr I ds{D j -D l ){T)g{T-s)S j {s) 



dr / ds ( Dj - t>i ) (r)g(r - s)Si(s) 



2 Jo 



1 .wAMk 
coth(— - — / 



4u, 



uj / dsg(s)(D j - Di)(s] 



(A.15) 
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General solution for arbitrary initial uncorrected states can be derived with the 
aid of a conjugate Wigner transform; we define with the initial system density matrix 
Pi(X = q + q' = q-q'), 

/oo 
dxe^ 2 p t (X,0- (A.16) 
-oo 

Tracing out the final variables leads to 

1 roo roo 

trp 0,0 / dx dt Pi (x,Z)tr J u ' l) = (A.17) 

2 J —oo J —oo 

exp [- jf* dr jf dsUr)P R (r - s)Co(s) + \ J* dr (S^r)^) + ^(t)6(t) ) 

•fx (jf ^^(r)(^(r) - A(r)),/ o * dr «?(t)( D,-(t) - Z>,(r) ) ) . (A.18) 

Functional differentiation with respect to j(r) ,/(r) gives the correlator such as 

(q(n)q{T2)) = ~-9(n - r 2 ) + / / ds^(ri - r)/3 fl (r - s)#(r 2 - s) 

2 jo Jo 

_ 9(Tl)9(T2) (^) M _ 9(Tl)g(T2) ^) M 

- (J7(r,)s(r 2 ) +j( Tl ) 9 (r 2 ))(|&) . (A.19) 

Here the suffix 00 is understood to mean p = , £ = . Needless to say, the derivatives 
are related to the averages of dynamical variables; 



i 



\ d ? J oo V d P J oo l^Woo 2 

(A.20) 

We assumed in deriving this formula that 

(fir ' (^L-°- (A - 21) 

The initial memory effect appears in two ways for the exact result of the correlator 
(g(Ti)g(r 2 )); first, via the initial state dependence, f p %, and secondly, via the explicit 
lower limit of the time integration, r = taken to be the initial time. The memory 
effect thus violates the time translation invariance under — > r, + <5. 

The coincident limit of the correlator is computed, using g(0) = 0, and the 
relation, 

rt rt roo 

/ dr ds g(t - t) (3 r (t - s) g(t - s) = du r + (uj) \h{uj ,t)\ 2 . (A.22) 
Jo Jo Jo 
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The result is eq.( |3.28|) and similar ones for the other quantities, eq.( |3.29|) , (|3.30|) in 
the text. 

We finally give some other examples of the conjugate Wigner function f p ^, 

1. thermal state of temperature T = l/(3 

UP , = exp[ -~ coth(^) ( V - + co o e ) ] , (A.23) 

2. moving packet of momentum po, with a spread Ap « ^Ja 

f p( (p , = exp[ ~ ( ^ + ai 2 ) + ipot) • (A.24) 
In these cases the correlators such as ([A.19Q can be worked out explicitly. 



Appendix B Renormalization of distribution func- 
tion 

One may use the equal time limit of the full propagator for the purpose of 
renormalization; in particular, 

( d d \ ~ 

hm — ) G(x ,y ;k) = Z , (B.l) 



+ \dx dy 0/ 



i ( d d N " 



hm - — - — G(x ,y ;k) = (H(k , t)) , (B.2) 
X0 ->y+ 4 \dx dyoj 

where Z is the wave function renormalization factor and TC(k,t) is the Hamiltonian 
for the momentum k mode. Using the expansion fl2.52|) , we find 

lim (■£--£-) G(x,y)= f e -*0HO (l +/(£,*)- /(£,*)) 

= 6%x-y), (B.3) 

in contradiction to eq. ( |B.1| ). Thus, the original expansion ( |2.52| ) must be modifid 
to allow a counter term in — iG(x ,y) . The term independent of / should thus be 
replaced as 



e 



iu k (x -y ) + g-twfctxo-w) ( 1 + / (fc ,t)) 



e iu> k {x - yo ) 1 1 _ + f(k tt j\ + e -i«>k(x -yo) + /(£,*)) . (B.4) 
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Another useful relation for the renormalization is the Fourier inversion formula, 

2 -i- (l + 2f(k,x )-2v(k,x )) 

/-> r°° dkn ~ -> 

d 3 xe^ x ( ¥ >(x,xoMO,y )) = -t / G(k ,k;x ) . (B.5) 

The correlator G(x ,y) = i{(f(x)ip(y)) in thermal equilibrium is what is called the 
real-time thermal Green's function in the literature |12|], and its Fourier transform 
G here is related to the analytic function ( [3.1 1| ) by 

- G(k , k) = 1 _ l e _ pkQ F(k + z0+ , k) + F(k - z0+ , k) . (B.6) 

This relation holds when ip particles are in thermal equilibrium with x- We shall 
first derive the renormalization condition assuming thermal equilibrium, and then 
extend its result to the non-equilibrium circumstance. 

Subtraction term for the distribution function consists of an infinite quantity 
and its associated finite term. In the spirit of the on-shell renormalization in field 
theory it is important to identify the pole contribution with possible infinities of the 
wave function and the mass correction included. The pole part of the propagator 
is extracted from the terms of F~ 1 (ko , k) to the quadratic order in ko', prior to the 
r enormalizat ion , 

iF polc (k ±iO + ,k) = , (B.7) 

kl — k 2 — M 2 (T) — n(fc , k) ± mr^(uj k , k) e(k Q ) 

u(.,k)=-v r d u' r -4^±=-v r ^'-vf , (b.8) 

J-oo OJ — UJ JO OJ i — 10* 

with M 2 (T) = M 2 + j^T 2 including the 0[X] temperature dependent mass. 

The mass and the wave function renormalization is done perturbatively for the 
proper self-energy U(ko , k). Expanding in powers of k 2 , — M 2 and k 2 , one identifies 
the renormalized temperature dependent mass as 

M 2 R (T) = M 2 (T) + n(M,0) , (B.9) 

and the wave function renormalization factor as 



''O / k =M 

Thus, defining the subtracted finite part 511 by 

5Tl(k , k) = Tl(k , k) - n(M , 0) + ( Z- 1 - 1 ) ( k 2 Q - M 2 ) , (B.ll) 
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one has 



; F pole (k ± *o + , k) = ^—^7 — A7 T- > (B.12) 

ko ~ (wf ) 2 ± wr_(^ fc , ft) e(ft ) 



a; 



ft 2 + Mr(T) + , ft) . (B.13) 



This F po } e (fco , k) is an optimal Gaussian approximation for the low energy dynamics, 
to the order 0[A 2 ]. One may regard SU(uJk , ft) as a finite energy shift due to the 



interaction with environment. When one uses UJ k = y ft 2 + M\ , for the energy at 
the pole, as is done in the on-shell Boltzmann equation, one has to compensate for 
the difference 



u?-u> k M m{ ? k ' k) , ^ = y^ + M 2 , (B.14) 

in the distribution function as a kind of finite energy renormalization. 

Relativistic covariance requires the equality of the temperature independent part; 



(B,5) 



with fcj a spatial component of ft. We shall check in Appendix D that the infinite 
part of this relation holds. The finite, temperature dependent part however needs 
not to satisfy the covariance relation, 

d 2 U(k ,k) d 2 U(k ,k) 

dk 2 = 9ft^' (R16) 

due to the presence of the preferred frame in which the temperature is uniquely 
defined. 

One may then compute (P-5| ) in the narrow width limit (or in the weak coupling 
limit). With v — 0, the renormalized pole term /^ le is given by 



2u k v ^ 2< 
We thus find that 

tR (. .R\ _ u k ~ _|_ u k 
2Zu)k Zlu 

from which to 0[X 2 ] 



1 + 2/ pole K)) = =-g (1 + 2f%M)) ■ (B.17) 



/pole(<) = + 1^7 UoleM , (B.18) 



SfUk) = /* ole (u;£) - / pole K 

5U{uj k ,k) \ fl \ 5U(u k , k) d/p ie 

~ 5Z + —2^) U + /polc J~^^^7' (R19) 
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with 5Z = Z — 1. Since the occupation number f(k) is defined in reference to the 
energy of the harmonic oscillator of mode k, (p\/ (2a;&) + Ukql/2 ) the finite part of 
the proper self-energy 511 (uk , k) appears in this formula as a change of the reference 
energy u k -> . 

In the low temperature limit the distribution function / po i e is very small, and one 
approximately has 

SUfivU-SZ+^ISM. (B.20) 
We extend this renormalization term to the case of the non-equilibrium state, to 

get 

1 ( ^ , 5U(uk,k) 

to be used at low temperatures. The renormalized distribution function is then 

8Z 5U(u)k , k) 
J 

The last term of this formula is related to the finite energy shift of the harmonic 
oscillator for the mode k. Its effect is absorbed by changing the reference energy 
when one defines the occupation number for this mode. In another word, this last 
term disappears by the replacement, u k — > , with uj} given by ( |B.13| ). Thus, if 
one modifies the occupation number according to 

2 R i 

fneAk)^(^ + ^qt)-^, (B.23) 

then the renormalized occupation number becomes 

A + Bk 2 

3 



SfUk) = -l-6Z+ 2 J 2 ' ' ) , (B.21) 



fUk) = f(k) + 5f vcn (k) « f(k) - — + V 4 J 2 ' ' ■ (B.22) 



/ren(fc) = f{k) , B = 8Z. (B.24) 

4kA 



Here A and B are infinite counter terms independent of the momentum k. 

Thie new definition appears very reasonable, because the ip system in isolation 
from the x environment is an ideal setting which has nothing to do with actual 
observation. The new reference energy uj} includes the interaction with the environ- 
ment, and in this sense is directly related to observation. We studied the problem 
more closely than this, by retaining the energy shift term 5U(uJk , k)/4w^.. We found 
a disease with this term; at very low temperatures the equilibrium distribution does 
not exist. This happens by having a dominantly negative term for the stationary 
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distribution / cq . Since this is clearly unacceptable, we shall use the definition ( |B.23| ) 
and its associated renormalization of the proper self-energy to define our occupation 
number. 



Appendix C Kinetic equation for unstable parti- 
cle decay 

We consider two-body decay of a boson ip described by a Lagrangian density, 

Ant = -^X 2 , (CI) 

where \l is a coupling constant of mass dimension. The environment particle \ is 
taken to make up a thermal bath of temperature T = 1/13. For simplicity, we assume 
that the mass of % particle vanishes. 

The fundamental quantity for the quantum kinetic equation of the unstable par- 
ticle decay is the spectral weight given by 

n 2 ( 1 1 1 — p-( 3uj + \ 

= < k °) { 2 *(M - k ) + Jk ln J ' (C - 2) 

u±=^(\k \±k), (C.3) 

where k = \k\. This spectral function corresponds to ip «-> XX f° r k > k and 
V^X ^ X f° r < ko < k, as depicted in Fig. 12. Note that we do not assume the on- 
shell kinematic condition k^ — k 2 = M 2 , thus the process (fx X becomes possible. 
For ko < the spectral weight is extended according to r x (— k) = —r x (k) . The 
Feynman- Vernon kernel 

a(x-y) = (^) 2 tT (T[x 2 (x)x 2 (y) P ,]) (C.4) 



given by eq. ( [2.31 ) contributes to the influence functional of this problem; 



F[<P > f'} = ex P[ - / dxdy(<p(x)-<p'(x)) (a(x-y)<p(y)-a*(x-y)<p'(y))}. 

Jx >yo 

(C.5) 

Neglect of higher powers of <p with the kernel of the type, (x 2 X 2 ' ' ' X 2 ) > defines the 
Hartree approximation in this decay model. 
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The Markovian kinetic equation can be worked out in the same way as in the 
annihilation-scattering problem in the text, and it is 



^^ = -r fc (f(k,t)-Uk)) 



(C.6) 



Uk) = f duj + |-) coth 



uj (3uj r x( u ' ^) 



2wjfc 2 (^ 2 -^)2 + ( 7rrx ( CJ; A;)) 2 

- ( T = contribution ) (C.7) 

*» T -&4 I + - — | ■ (Gf 



Here we ignored a minor correction, the temperature dependence of mass, and took 
an advantage of the weak coupling limit in writing the last form of this equation. 
The rate Tk = tt r x (u>k , k) / u>k is the decay rate of unstable ip particle with time 
dilatation effect included. 

As demonstrated in ||, the equilibrium distribution function can be interpreted 
by a Gibbs formula, 

where H tot is the total Hamiltonian including interaction between the ip system and 
the x environment. 

— * 

Let us work out the stationary distribution f eq (k) in more detail. The narrow 
width approximation to the energy integral (|C.7|) gives the the usual Planck dis- 
tribution, and the rest of contribution at low temperatures is approximately given 

as 

1 r A,, r x(u,k) 



oji Jo e 13 ^ — 1 



k 

2 roo 1 / 1 1 1 - p-0(w+fc)/2 \ 

5-^/ da; -= - 0(w - A;) + — In _, . (C.ll) 

16tt 2 Jo - 1 V 2 V ; 1 - e-^-m J K J 

The last u integral can be worked out analytically both in the limit of k <C T and 
k ^> T. For fc<T the dominant part of the integral extends both from the region 
uj < k and to uj > k, while for k 3> T it comes only from uj < k, giving 

2T C(2)T 2 1 

e */2T_l' fc e */2T_l' 
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respectively. Here ((2) = ^- is the Riemann's zeta function. Our interpolating 
formula is a smooth match of these two limiting functions, 

Jcq{ ' 16tt 2 u>\ 2k + C(2) T e k / 2T - 1 ' 1 ' 

We numerically compare this interpolating formula with the result of exact in- 
tegration at T = M/20 in Fig.13, where we plot the quantity k 2 f eq (k) / (2n 2 ) using 
( |C.7|) . The agreement of our approximate formula f eq (k) (|C.13|) and the exact nu- 
merical integration is not excellent, but is adequate for evaluation of the total num- 
ber density. Thus, at low temperatures the stationary distribution for the unstable 
particle is approximately given by our 5f cq (k), eq.( |C13|) . 

Integration over the phase space gives the number density at low temperatures 
of T <C M (parent mass); 

Sn eq = At r dk k 2 5f cq (k) = C( ?/f, cT\ (C.14) 
q 2vr 2 Jo J qV ; 2vr 4 M 3 K J 



oo <t»2 



dx ^-—^^^0.27. (C.15) 

Ax + C(2) e x -l K J 



Thus, 

MT T 
n eq (T) * (^) 3 / 2 e- A ^ + 0.23r(^) 2 T 2 , (C.16) 

where V = fi 2 / (32ttM) is the decay rate for if — > XX m t ne V 9 res t frame. See Fig. 14 
for the number density. Emergence of the temperature power term has been found 
recently fl5f , || , and its relevance to the GUT baryogenesis has been discussed in [|7| . 



Appendix D Computation of off-shell distribu- 
tion function 

We separate both 6f eq into three pieces and write 

6f eq + 6f mn = h + f 2 + / 3 - A \ B f , (D.l) 

4a; fc 

= -L r^ ^l^| , (1=1,2), (D.2) 
2^ J-oo (a; - cj fc ) 2 + Iy4 

r f,, A-2 / d 3k ' r x {\u-u k ,\,k-k') p 

ri( ^' fc) - 2 i (2vr)3 2 ^ e^l—vl-1 /( *°' (D ' 3) 
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^^= 2 / (x. rx( '^'l-l +y) (/(g)+1) - (D4) 

We shall give the explicit form of fz{k) later, which has the temperature dependence 
only via the two-body spectral function r x . For each piece we shall give an integrated 
number density to show their respective importance at low temperatures. 

The first piece fi has the main contribution from the region, \uj — uj k >\ < T, due to 
the exponential suppression outside this region. One may use the expansion formula, 

r x (\u-u k ,\,h-k') ( dr x (x,k-k') \ 

\ / x=0 

which gives the corresponding spectral function, 

This contributes to the stationary distribution function off the mass shell; 

A(*>3££- r**-f& r ^f^-/^\ (D.7) 
327T 4 ku k Jo <jj k i J-oo [uj — oj k y + 

1 _ e -p(k+k')/2 

D(u,k,k') = 9(\k - k'\ - \u - u k >\) In - _ e _ /3|fc _ fc/|/2 

1 _ e -/3(fc+fc')/2 

+ - u*| - |* - In j^z^z^ • (D.8) 

Since the k! integral is dominated in the region around k of width of order T, the 
integral roughly gives 



x /(*)■ ( D - 9 ) 



W ^°W X ^H? n ^ (D - 10) 



Integrated over momenta, this gives the number density 

A 2 T 2 

J k 

with n v the tp number density. One may interpret this result, by saying that a ip 
fraction of order \ 2 T 2 /M 2 is created by scattering off the mass shell. 

Since the magnitude of the first term n{ is important to estimate the relic abun- 
dance, we performed direct numerical computation for this term without resorting 
to the above approximation. The number density due to this term is 

f "f^f™ ,mcM ' (D - U) 
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cm = r ix r^" *» f da ,^+")+^-"+")-^w ;(D . 12) 

X 1 _ e -|w+a;|/2 

A(u) = In : m , (D.13) 

in an energy unit measured by the temperature T. The narrow width limit r —>■ can 
trivially be taken in this form of the integral. We numerically confirmed that C(k') 
is alomost in proportion to k' such that the last k' integral approximately gives the 
tp number density n^, and moreover the constant factor is roughly 5.5 x (T/M) 0-35 
within an accuracy of 10 %, if the temprature range is M/100 < T < M/10. This 
gives 

5n{ « 5.5 x ^ nv « .89 x 1(T 5 A 2 (^) 1J » »„ . (D.14) 

This result has a slightly slower decrease in temperature than the analytic estimate 
of oc T 2 n ¥ . It nevertheless gives a subdominant contribution for n eq at low temper- 
atures. 

The second piece f 2 is dominant in the region around u ~ — cu^ ; 

, /A 2 /• <i 3 /c' 1 f°° , rJuo,k + k') \ , 

Uk) « — / - — / ck> *V ; ( /(A;') + 1 ) D. 15 

= fi(k) + f 2 (k) . (D.16) 

We separated the second contribution into a sum of /—dependent and /—independent 
integrals. The last u integral here is explicitly done in Appendix C in a related model, 
giving with q = \k + k'\ 

r^g) 4 C (2)T 2 1_ 

This gives, for the /—dependent part, 



16tt 4 7o uu {uj k + u k ,) 2 J\ k -k'\ 2q + C(2)T e<?/ 2T - 1 ' 

(D.18) 

Considering that the q integral is suppressed at q ^> T, one obtains for the integrated 
number density, 

for the / dependent part. Here c ~ 0.27 from eq. (|C.15| ) of Appendix C. This is 
smaller by a factor of order (T/M) 2 than the first piece n[. 
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On the other hand, the / independent part gives for the distribution function 
and the number density, 



n. 



167T 4 ku k Jo 2q + C(2)T e^ 2T - 1 \u k + u k - q u k + u k+q )' 

(D.20) 

The total number density thus derived is roughly of order A 2 T 4 /M. 

Physical processes that contribute to the second piece /° are predominantly in- 
verse annihilation XX ~^ <Wj an d 1 to 3 process, x ~^ XW-> which gives a small 
fraction of the number density. 

The third piece fz{k) has the following spectral function, 

r 3 (uj,k)=2 J ^ 2ix y 2uJki ~ u) r x(\^ ~^k'\,k- k') f(k') 

+ 6(-Lu k , -Lu)r x (\u + uj k/ \,k + k>)(l + f(k>))) , (D.23) 



\2 / 2 1 - e'^H+i)/ 2 \ 

One can separate the temperature independent (f^) and dependent (/^) pieces 
for the third contribution by the presence of the (3 factor; prior to the subtraction 
of the counter term, 



128tt 4 koj k J ~'" oj k > J-oo (uj - uj k ) 2 + T|/4 

rk+k' 

■ ^ dqq {Gi{u ,u k > ,q) - Fi(u k > -u k ,q) f{k')) , (D.25) 

Gi(u , u k , , q) = Fi(- u - w, q)(l + f(k') ) + F,(- cu + w w , q) f(k') , (D.26) 

2T 1 - e ~^ x \ +q)/2 

F (x , q) = 9(x - q) , F T (x , q) = — In f — -^-^ . (D.27) 



There are divergent terms of the form, 

A + Bk 2 



*4 



(D.28) 



which are cancelled by the mass and the wave function counter terms in the proper 

— * 

self-energy II (a;, A;). To see this, let us focus on the term both independent of the 
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temperature and of the distribution function /, which can be worked out by the u 
integration explicitly; 

/f^fcj^-i- [ A dk> 
Js v ' 128vr 4 ku k Jo 

■— ( k + k'-\k-k'\- {u w + u> h ) In Uk ' +Uk + k + f ) (D.29) 
cd k ' V + tUfc + | A; — ft' | y 

where A is the momentum cutoff. By expanding around the momentum k = 0, we 
may identify 

A = - M / tZJf 7 — , (D.30) 

16vr 4 7 (w fc / + fc' + M) v ; 

A A 2 r A fr' 2 
5 = ^--^— / tZA/ — -. (D.31) 

2M 2 32vr 4 Jo u k , {uj k , + k' + M) 2 v ; 

The logarithmic infinity, the second term of eq. (p.31 | ) , is related to the two- loop 
infinity of the self-energy diagram of Fig. 15. By a straightforward calculation this 
diagram gives a wave function renormalization factor, 

\2 

5Z= In A 2 . (D.32) 

256vr 4 v ; 



One sees that this corresponds to the logarithmic infinity of eq. ( p. 31]) , confirming 
the relativistic covariance of the wave function renormalization. 

Removal of infinities thus works with renormalization. We are not much inter- 
ested in the remaining, temperature independent finite part, since this is absorbed 
by a new definition of the vacuum in interacting field theory. 

(T) 

We thus next turn to the finite, temperature dependent part o/ 3 . For some 
part of this integral it is not difficult to count the mass M dependence in the low 
temperature limit. Thus, the integrated number density from this part is either of 
the form, A 2 (T '/ 'M) 4 n v or A 2 (T/M) 4 T 3 . These are smaller by a positive power of 

than 5rii and 5ri2 we already considered. The remaining piece has the spectrum 
9((jJk> — u})r!jp(\u — uy \ ,k — k') with r*p the temperature dependent two-body 
spectrum. This contribution can be combined with f\(k) of eq.( p.7] ), to give the 
following /—dependent distribution, 

_ X 2 T r°° k' rk+k' roc 

//(*) = ^riT- dk ' — f{k,) I , dq / du 

D47T 4 kUk JO UJk' J\k-k'\ J~oo 

1 (l 1 _ e -/3(|o;-uv|+5)/2 \ 

ln \ -«\\ l ,-,„..\-n\/7 - (w -> . (D.33) 



u) k ) 2 + r 2 k /4 \ e ^-"v)-l 1 - e-/ 3 ! 1^-^1-9 1/2 
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We have numerically computed this k' integrand, which turns out well in proportion 
to k' 2 in the non-relativistic limit. Thus, the integrated number density is of order, 
A 2 (|j) 2 , of the same order as n[ of eq. ( p,10| ). 

In summary, the dominant term of the off-shell contribution is 

Sf eq (k) « f f (k) + f 2 (k) , (D.34) 



where // is given by eq.( p.33|) or its numerically better alternative, and /° by 
eq.( p.2Uj ). The use of the more precise numerical result for the /—dependent term ff 
such as ( p.14 ) is not necessary, since it is shown in Section 4 that the /—dependent 



term f® dominates in the evolution equation for the number density, thus 5f eq ~ f® ■ 
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Figure caption 



Fig.l 

Tadpole diagrams. 

Fig.2 

Two-body spectral function r x (u , k) for two massless XX state. Two choices of 
the momentum k relative to the temperature are shown. 

Fig.3 

Four different processes contributing to the Boltzmann equation. The solid lines 
are for the p> particle, while the dotted are for the \ particle. 

Fig.4 

Singularity structure of the self-energy function F(z). The dotted cross gives 
poles in the second Riemann sheet, while the wavy line is the branch-cut singularity. 

Fig.5 

Functions of the integrand in the energy integral, eq.( p.57|) in the text. The 



dotted line is the Planck distribution of T = 0.01 M, while the dashed is the rest of 
the Breit-Wigner form, with the solid line giving product of these two. Parameter 
values taken are k = 0.01 M ,k' = 0.005 M , A = 0.01 . 

Fig.6 

Time evolution of the yield Y for <p mass of 10 GeV, and indicated couplings. 
For a comparison the on-shell Boltzmann result is also shown. 

Fig. 7 Time evolution of the yield Y for the A = 10 -3 and a few p> mass values. 

Fig.8 

Comparison of the S-wave and P-wave annihilation in the time evolution. 
Fig.9 

The parameter region for the off-shell and the on-shell dominance in the S- 
wave boson pair annihilation model. The contour of the closure density, p c = 



58 



10 4 eV cm 3 , is also shown. The parameter relations are 5 = 5.4 x 10 6 A 2 ,5/rj = 
1.2 x 10~ 22 M/GeV . 

Fig.10 

Contour lines for the present mass density in the unit of p c = 10 4 eV cm~ 3 (S-wave 
model). 

Fig.ll 

Similar contour lines to Fig.10 for the P-wave annihilaation. 
Fig.12 

Diagram for two-body spectral function r x . 
Fig.13 

Equilibrium distribution function for the unstable particle decay at T = 0.05 M. 
The exact computation (solid) is compared with the approximate one (dotted), and 
the Planck formula (dashed). In the inset result is given in linear scale. 

Fig. 14 

Equilibrium number density; exact result (solid) is compared with the on-shell 
result (dotted). 

Fig.15 

Self-energy diagram. 
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